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Executive  Summary 


The  research  described  in  this  report  was  conducted  in  fulfillment  of  Project  MM- 1592, 
“Non-traditional  physics-based  inverse  approaches  for  determining  a  buried  object’s  location,” 
submitted  to  the  Strategic  Environmental  Research  and  Development  Program’s  (SERDP) 
Exploratory  Development  Program  (SEED)  in  response  to  Statement  of  Need  MMSEED-07-02, 
“Advanced  Technologies  for  Detection,  Discrimination,  and  Remediation  of  Munitions  and 
Explosives  of  Concern  (MEC):  UXO  Technology.” 

The  main  focus  of  this  research  was  to  explore  and  develop  new  electromagnetic  inverse¬ 
scattering  approaches  for  estimating  the  locations  and  orientations  of  buried  objects  that  would 
not  involve  the  solution  of  an  ill-posed  nonlinear  minimization  problem  and  that  would  be 
sufficiently  robust  and  efficient  to  be  used  in  real  field  UXO  discrimination.  We  concentrated  on 
the  fundamental  aspects  and  potential  practical  applicability  to  the  UXO  problem  of  three  novel 
physics-based  approaches:  (1)  pole-series  expansions,  (2)  energy  gradients,  and  (3)  use  of  the 
physical  properties  of  left-handed  media  (LHM).  All  three  methods  assume  that  targets’ 
responses  can  be  approximately  reproduced  by  means  of  a  set  of  magnetic  dipoles — the  scattered 
field  singularities — distributed  in  particular  points  inside  the  objects.  The  algorithms  are 
implemented  using  a  numerical  technique  called  the  Normalized  Surface  Magnetic  Source 
(NSMS)  model.  We  studied  the  accuracy  with  which  the  methods  can  estimate  an  object’s 
location,  orientation,  and  equivalent  magnetic  polarizability,  their  robustness  with  respect  to 
noise,  their  computational  speed,  and  their  requirements  with  regard  to  data  quality  and  quantity. 
We  ultimately  wanted  to  determine  which  of  these  three  techniques  would  be  practical  and 
reliable  for  use  in  real  UXO  production  runs. 

As  part  of  this  investigation  we  studied  several  existing  data  sets  collected  by  established 
instruments  such  as  the  Geophex  GEM-3  and  the  Geonics  EM-63  and  by  next-generation  sensors 
like  the  NRL  discrimination  array,  the  MPV-TD  and  the  GEM-3D+.  The  GEM-3,  GEM-3D+,  and 
TD  MPV  data  were  collected  at  the  ERDC-CRREL  test  site  in  Hanover,  NH,  by  CRREL 
personnel;  the  EM-63  data  were  collected  at  the  ERDC  test  stand  in  Vicksburg,  MS,  by  operators 
from  Sky  Research  and  the  University  of  British  Columbia;  and  the  NRL  discrimination  array 
data  were  collected  by  NRL  researchers.  The  new  proposed  algorithms  were  applied  to  these  data 
sets  in  order  to  estimate  locations,  orientations,  and  equivalent  magnetic  polarizability  tensors  of 
buried  objects;  in  addition  we  conducted  discrimination  studies  using  the  NSMS  model  so  as  to 
demonstrate  the  applicability  of  these  techniques. 
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Chapter  1 
Introduction 


The  research  described  in  this  report  was  conducted  in  fulfillment  of  Project  MM- 1592, 
“Non-traditional  physics-based  inverse  approaches  for  determining  buried  objects’  location,” 
submitted  to  the  Strategic  Environmental  Research  and  Development  Program’s  (SERDP) 
Exploratory  Development  Program  (SEED)  in  response  to  Statement  of  Need  MMSEED-07-02, 
“Advanced  Technologies  for  Detection,  Discrimination,  and  Remediation  of  Munitions  and 
Explosives  of  Concern  (MEC):  UXO  Technology.” 

Several  data-processing  techniques  for  geophysical  survey  data  have  been  developed  for 
discriminating  between  unexploded  ordnance  (UXO)  and  non-UXO  items.  Typically  the  first  step 
of  these  methods  is  the  recovery  of  a  set  of  parameters  that  specify  a  physics-based  model 
representing  the  object  under  interrogation.  For  example,  in  electromagnetic  induction  (EMI) 
sensing,  the  recovered  parameters  consist  of  the  object’s  location  and  spatial  orientation  in 
addition  to  “intrinsic”  parameters  such  as  the  polarizability  tensor  (along  with  some 
parameterization  of  its  time-decay  curve)  in  dipole  models  (Pasion  and  Oldenburg,  2001)  or  the 
amplitudes  of  responding  magnetic  sources  in  the  NSMS  and  SEA  models  (Shubitidze  et  al , 
2005;  Sun  et  al .,  2005;  Shubitidze  et  al .,  2007).  EMI  responses  depend  nonlinearly  on  the 
subsurface  object’s  location  and  orientation,  and  inverting  for  these  parameters  quickly  and 
accurately  is  a  crucial  part  of  all  the  inversion  techniques  currently  in  use  by  the  UXO 
community. 

Most  EMI  sensors  are  composed  of  separate  transmitting  and  receiving  coils.  When  the 
operator  activates  the  sensor,  a  current  runs  through  the  transmitter  coils,  which  results  in  the 
establishment  of  a  (“primary”  or  “principal”)  magnetic  field  in  the  surrounding  space.  By 
Faraday’s  law,  this  time-varying  magnetic  field  induces  eddy  currents  in  highly  conducting 
bodies;  ferromagnetic  bodies  also  have  their  magnetization  affected  by  the  impinging  field.  These 
currents  and  magnetization  in  turn  generate  a  (“secondary”  or  “scattered”)  magnetic  field  that 
also  varies  with  time  and  induces  measurable  currents  in  the  receiving  coils.  At  the  end,  the 
electromagnetic  data  are  inverted  using  different  forward  models.  The  procedure  for  estimating 
the  location,  orientation,  and  electromagnetic  parameters  of  a  buried  object  (strung  together  in  a 
“model  vector”  v)  is  carried  out  by  defining  an  objective  function  that  quantifies  the  goodness-of- 
fit  between  the  measured  data  and  the  predictions  of  the  forward  model.  Routinely,  a  least- 
squares  (LS)  approach  is  taken  to  recover  v:  formally,  if  dobs  is  the  vector  of  the  measured 
scattered  field  and  F(v)  is  the  solution  to  the  forward  problem,  the  least-squares  criterion  assumes 
the  form 


minimize  4>(v)  =  dobs  -  F ( v) 


(1) 


A  simple  way  to  determine  the  model  vector  v  is  to  use  the  Gauss-Newton  method,  which  starts 
with  an  initial  guess  v0  and  updates  it  iteratively  through 


(2) 
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where  k  denotes  the  iteration  number  and  is  a  perturbation  direction;  we  solve  for  the  that 
minimizes  </>.  This  approach  is  computationally  intensive  because  it  requires  massive  repetition  of 
the  forward-problem  solution.  In  particular,  the  determination  of  the  nonlinear  elements  of  v — the 
location  and  orientation  of  the  object — is  nontrivial  and  time-consuming. 

The  main  motivation  of  this  project  was  to  develop  new  non-traditional  physics-based 
inverse  scattering  approaches  that  would  enable  us  to  infer  the  location  of  a  buried  object  from 
measured  data  without  having  to  solve  a  nonlinear  ill-posed  EMI  inverse  problem.  Specifically, 
we  proposed  to  investigate  three  physics-based  methods  involving  (1)  energy  gradients,  (2)  pole- 
series  expansions,  and  (3)  the  physical  properties  of  left-handed  media  (LHM).  All  three 
approaches  assume  that  targets  exhibit  a  dipolar  response.  Secondary  magnetic  fields  in  the  EMI 
regime  are  due  to  eddy  currents  or  magnetic  dipoles  distributed  nonuniformly  inside  the 
scatterers.  There  are  some  particular  points,  the  “scattered  field  singularities”  (SFS),  where  most 
of  these  sources  are  concentrated.  Recent  studies  show  that  under  certain  conditions  the  entire 
scatterer  can  be  replaced  with  several  responding  elementary  sources  placed  at  the  SFS 
(Kyurchkan  et  al. ,  1996;  Zaridze  et  al. ,  2002;  Bliznyuk  et  al. ,  2005).  The  mathematical  and 
physical  properties  of  SFS  and  their  application  to  electromagnetic  scattering  problems  are  very 
well  documented  and  studied,  and  constitute  the  discipline  known  in  the  literature  as  “Catastrophe 
Theory”  (Arnold,  1986).  Our  objective  in  this  project  was  to  utilize  the  three  approaches 
mentioned  above  to  find  fast  and  effective  ways  to  determine  the  locations  of  these  singularities. 

In  the  energy-gradient  approach  we  found  that  it  is  possible  to  determine  the  location  and 
orientation  of  a  buried  object  if  we  have  access  to  (1)  the  magnetic  field  vector  H,  (2)  the  vector 
potential  A,  and  (3)  the  scalar  magnetic  potential  y/  at  a  single  point  in  space.  Since  of  these 
quantities  only  H,  or  one  of  its  components,  is  available,  as  part  of  this  project  we  developed  a 
numerical  procedure  based  on  the  2D  normalized  surface  magnetic  source  (NSMS)  model  that 
reconstructs  A  and  ^from  H. 

The  NSMS  model  can  be  considered  as  a  generalized  surface  dipole  model  and  indeed 
reduces  to  the  usual  dipole  model  in  a  special  limiting  case.  The  full  3D  NSMS  approach  models 
an  object’s  response  to  a  primary  field  by  distributing  a  set  of  equivalent  elementary  magnetic 
sources  over  an  auxiliary  surface  surrounding  the  target  (Shubitidze  et  al .,  2007);  the  field 
radiated  by  these  elementary  sources  satisfies  analytically  the  governing  EMI  equations  outside 
the  scatterer.  The  surface  integral  of  the  source  amplitudes  can  be  used  as  a  discriminator; 
moreover,  the  method  can  be  used  to  generate  a  synthetic  scattered  magnetic  field  H  and  its 
corresponding  potentials  A  and  y/.  The  NSMS  has  proved  to  be  successful  in  single-  and 
multiple-target  discrimination  studies,  but,  lacking  a  method  to  rapidly  and  independently 
determine  the  location  and  orientation  of  the  target  (and  hence  of  the  surrounding  surface),  it  has 
to  be  used  within  a  time-consuming  and  computationally  expensive  optimization  procedure.  The 
2D  NSMS  method,  on  the  other  hand,  replaces  the  surface  around  the  scatterer  with  a  flat  plane  of 
dipoles  at  a  (known)  location  intermediate  between  the  instrument  and  the  target.  Given  high- 
spatial-coverage  geophysical  data  taken  over  a  buried  object  it  is  possible  to  calculate  the 
amplitudes  of  the  responding  magnetic  sources  on  the  intermediate  plane  by  solving  a  linear 
system  of  equations.  The  amplitudes  can  then  be  used  to  estimate  H,  A,  and  ^at  any  point  on  or 
above  the  measurement  surface,  and  these  in  turn  can  be  combined  analytically  to  solve  for  the 

relative  location  R  and  the  polarizability  M  of  the  hypothetical  dipole. 

Chapter  2  outlines  our  research  into  this  method,  to  which  we  will  refer  from  now  on  as 
the  HAP  approach.  As  we  shall  see,  the  HAP  method  provides  fast  and  accurate  estimates  of  the 
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location,  orientation,  and  induced  magnetic  polarizability  of  a  visually  obscured  object.  We  tested 
the  procedure  using  both  synthetic  and  measured  data,  both  in  the  time  domain  (with  the  EM-63, 
MPV-TD,  and  NRL  array  sensors)  and  in  the  frequency  domain  (using  the  GEM-3D+),  and 
successfully  combined  it  with  the  3D  NSMS  method  in  MPV-TD  and  GEM-3D+  blind  tests.  Its 
performance  and  its  tolerance  to  noise  make  this  approach  our  definitive  candidate  for  real-field 
applications.  We  also  introduce  an  equivalent  but  reduced  version  of  the  HAP  that  employs  only 
H  and  y/. 

The  pole-series-expansion  approach  involves  determining  a  spatial  distribution  of 
scattered  field  singularities  which  can  be  related  to  a  target’s  properties  (/.<?.,  its  shape  and  size).  It 
uses  scattered-field  information  at  a  given  set  of  points,  again  estimated  via  a  2D  NSMS 
approximation,  to  quickly  locate  the  origin  of  this  scattered  field  at  any  spatial  resolution.  A 
similar  technique  developed  for  this  project  is  the  pseudospectral  finite-difference  (PSFD) 
method,  which  is  based  on  spatial  finite-difference  approximations.  Chapter  3  outlines  our 
research  into  these  modeling  techniques.  We  find  that  the  methods  are  fast  but  very  sensitive  to 
noise,  which  detracts  from  their  potential  usefulness  in  the  field. 

The  last  approach  we  developed  uses  left-handed  media  (LHM),  and  in  particular  their 
EM  field-refocusing  property,  as  a  tool  for  estimating  the  location  of  a  buried  object.  In  the  low- 
frequency  EMI  regime  one  first  reconstructs  a  Gaussian  beam  by  combining  a  monostatic  EMI 
spatial  distribution  with  an  empirically  assigned  wave-phase  part;  after  that  one  solves  a  full  EM 
scattering  problem  in  order  to  map  the  field  distribution  back  to  its  origin.  Chapter  4  outlines  our 
research  into  LHM  in  order  to  assess  their  potential  usefulness  in  the  context  of  the  UXO 
problem.  Our  studies  show  that  in  the  low-frequency  regime  the  LHM  refocusing  approach  is 
impractical  because  the  necessary  wavelengths  are  much  larger  than  the  size  and  depth  of  any 
potentially  interesting  target;  however,  the  research  reported  here  could  be  useful  for  the  EM- 
wave  community. 

Finally,  Chapter  5  summarizes  this  report  and  suggests  future  research  directions. 
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Chapter  2 


The  HAP  approach  for  estimating  the  location, 
orientation,  and  magnetic  polarizability 
of  a  buried  object 


2.1.  Introduction 

In  this  chapter  we  present  a  new,  physics-based  expression  for  determining  the  location, 
orientation,  and  magnetic  polarizability  of  a  buried  object.  The  approach  assumes  that  the  target 
exhibits  a  dipolar  response;  it  requires  three  global  values — the  magnetic  field  vector  H,  the  vector 
potential  A,  and  the  scalar  magnetic  potential  y/ — to  be  evaluated  at  a  single  location  in  space.  Of 
these  values  only  the  scattered  magnetic  field  H  is  measurable.  Therefore,  in  order  to  estimate  the 
scattered  magnetic  scalar  and  vector  potentials  from  data,  we  employ  a  numerical  technique  called 
the  normalized  surface  magnetic  source  (NSMS)  method.  In  the  original  NSMS  model,  the 
scattered  magnetic  field  outside  an  object  is  reproduced  mathematically  by  equivalent  magnetic 
charges  distributed  on  a  3D  closed  surface.  Here  we  use  a  2D  implementation  of  the  NSMS  that 
uses  elementary  magnetic  dipoles  (instead  of  magnetic  charges)  distributed  on  a  planar  surface 
placed  under  the  measurement  grid.  These  sources  are  used  to  estimate  the  scattered  magnetic 
field’s  vector  potential  A  and  scalar  magnetic  potential  y/  without  a  priori  knowledge  of  the 
object’s  location  and  orientation.  The  amplitudes  of  the  NSMS  are  determined  by  matching  the 
measured  magnetic  field  with  the  NSMS  modeled  field.  Once  the  NSMS  amplitudes  are 
determined,  H,  A,  and  y/  are  simulated  on  or  above  the  measurement  grid.  After  reviewing  the 
theoretical  basis  and  the  practical  implementation  of  the  new  technique,  we  demonstrate  its 
applicability  to  UXO  discrimination  by  using  synthetic  and  measured  data.  We  then  present  several 
numerical  and  experimental  tests  with  actual  UXO  interrogated  by  the  EM-63,  MPV-TD,  NRL  EMI 
array  and  GEM-3 D+  sensors. 


2.2.  Estimating  a  buried  object’s  location,  orientation  and  magnetic 
polarization  from  EMI  data 

In  the  low-frequency  EMI  regime  considered  here,  the  primary  magnetic  field  penetrates  a 
metallic  object  to  some  (frequency-dependent)  degree  and  induces  eddy  currents/magnetic  dipoles 
within  it.  These  induced  currents  and  dipoles  then  produce  a  secondary  or  scattered  field  outside  the 
object  that  is  measured  by  a  receiver.  Using  the  magnetic  dipole  approximation  we  can  write  the 
magnetic  field  H  and  the  scalar  and  vector  potentials  y/  and  A  as 


JkR 


H  =  - 


AnRi 


3R(R  m) 

R2 


-m 


(l  -  jkR)  -  &2  (R  x  (R  x  m)) 


(2.1) 
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Figure  2.1.  Geometry  of  the  problem.  The  target  is  assumed  to  be  a  magnetic  dipole  located 
at  rd. 


¥  = 


(R  m) 

4  nR3 


(1  -jkR)e», 


(2.2) 


a  =  M„ 


(mxR) 
4  nR3 


(1  -jkR)eikR, 


(2.3) 


where  k  is  the  wave  number  in  the  surrounding  medium,  R  =  r  -  rd  ,  and  the  vectors  r  and  xd  are  the 
observation  point  and  the  location  of  the  source,  which  we  assume  to  be  a  dipole  (Stratton,  1941, 
Chapter  8)  (see  Figure  2.1).  From  (2.1)  note  that,  in  the  electromagnetic  wave  regime,  the  magnetic 
field  due  to  a  magnetic  dipole  has  terms  that  decay  as  R~\  R  2,  and  R 3.  The  range  kR  » 1  is  referred 
to  as  the  far  zone,  and  fields  in  this  range  are  referred  to  as  being  in  the  far  field.  Similarly,  fields  in 
the  near  zone  (with  kR  « 1 )  are  referred  to  as  being  in  the  near  field,  while  the  zone  kR  «  1  is  called 
the  intermediate  zone.  Typically,  UXO  detection  and  discrimination  are  conducted  in  the  near  and 
intermediate  zones  during  both  land-based  and  underwater  UXO  cleanup.  Additionally,  in  the  EMI 
regime,  displacement  currents  are  considered  irrelevant,  which  means  that  the  contribution  of  the  k 2 
term  in  (2.1)  can  be  set  to  be  zero.  Under  these  assumptions  we  can  take  the  dot  product  of  (2.1) 
and  R  and  use  (2.2)  to  show  that 
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R3 
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G(R,k)  =  2i// , 
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where  G(R,k )  =  (l  -  jkR)ejkR  I  An  . 

Similarly,  we  can  take  the  cross  product  of  (2.1)  and  R  and  use  (2.3)  to  obtain 


HxR  =  G(R,k)— 
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Now,  taking  the  cross  product  of  H  and  (2.5)  gives 


Hx 


-A 


=Hx[HxR]  =H(H  R)-R  |H|2  =2H  ^-R  |H 


(2.5) 


(2.6) 


and  from  this  we  can  solve  for  R: 
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(2.7) 


2H  ^  +  [HxA /jUo' 


Thus  if  we  assume  that  the  source  is  a  single  dipole  we  can  express  its  location  vector  R  in 
terms  of  only  three  global  values.  Moreover,  this  expression  is  independent  of  the  frequency,  which 
means  that  (2.7)  is  valid  for  both  free  space  and  conducting  media  such  as  water  as  long  as  the 
MQS  assumption  holds;  i.e.,  as  long  as  we  can  neglect  displacement  currents.  Note  that  R  is 
determined  as  a  ratio  between  the  three  global  values.  This  makes  the  expression  in  (2.7)  partially 
tolerant  to  noise  due  to  scaling  arguments,  since  A  and  ^depend  on  the  H  field  ( cf.  (2.4)  and  (2.5)). 
The  magnetic  dipole  moment  can  be  found  by  taking  the  cross  product  of  R  and  (2.3)  and  using 
(2.2): 

m  =  ^f^(Rr-[A,ftxR]),  (2.8) 

where  R  is  determined  from  (2.7)  (see  Figure  2.1).  Note  that,  for  an  isolated  dipolar  source,  only  a 
single  instance  of  H,  A,  and  y/  is  needed  to  determine  R  and  m  uniquely. 


2.3.  A  reduced  HAP  formulation  for  estimating  a  buried  object’s 
location  and  orientation 

The  analytic  expressions  (2.7)  and  (2.8)  derived  in  the  previous  section  require  the 
magnetic  field  vector  H,  the  vector  potential  A,  and  the  scalar  magnetic  potential  yr  at  a  single 
location  in  space;  the  method  uses  seven  global  values  to  estimate  the  object’s  location,  which,  as 
we  can  see  from  (2.7),  consists  of  only  three  unknown  parameters.  We  now  present  a  reduced  HAP 
formulation  that  uses  only  the  magnetic  field  and  the  scalar  potential,  thus  decreasing  the  number  of 
required  global  values  and  reducing  computational  requirements.  Recalling  that  R  =  r  -  rd,  where  r 
and  rd  are  respectively  the  observation  point  and  the  dipole  location  in  the  global  coordinate  system 
(see  Figure  2.1),  we  can  write  (2.4)  as 


H  •  (r  —  rd )  =  2y/ 

(2.9) 

H  •  rd  =  -2y/  +  H  •  r  . 

(2.10) 

Since  r,  H  and  y/  are  known  at  every  observation  point,  the  last  equation  can  be  rewritten  in  matrix 
form  as  follows: 
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Under  the  low-frequency  EMI  assumption  k2  — »  0 ,  equation  (2.1)  can  be  rewritten  as 


(2.12) 
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which  lets  us  find  the  induced  magnetic  dipole  moment  via 


<2-i3) 

once  we  have  determined  R  from  (2.1 1).  Thus  the  location  rd  and  the  dipole  moment  m  of  a  single 
isolated  dipolar  source  can  be  uniquely  determined  from  H  and  y/. 


2.4.  Predicting  the  vector  and  scalar  potentials  from  EMI  data 

As  we  have  stated  before,  of  the  three  global  values  (H,  A,  y/)  needed  to  find  location, 
orientation,  and  magnetic  polarizability  we  have  access  only  to  the  H  field  (either  one  or  all  three 
components).  Therefore  we  need  to  find  a  way  to  obtain  the  potentials  A  and  ^from  the  measured 
magnetic  field  H  to  be  able  to  use  our  method. 

Recently,  new  EMI  sensors  have  been  developed  that  can  operate  in  both  monostatic  and 
multistatic  data-collection  regimes.  Monostatic  sensors,  such  as  the  Geophex  GEM-3  frequency- 
domain  instrument  (Won  et  al .,  1997)  and  the  Geonics  EM-61  and  EM-63  time-domain  instruments 
(McNeill  and  Bosnar,  1996),  have  collocated  transmitter  and  receiver  coils.  Multistatic  sensors  like 
the  Man-Portable  Vector  (MPV)  time-domain  instrument  from  CRREL  (Barrowes  et  al ,  2006),  the 
Berkeley  UXO  Discriminator  (BUD)  system  (Gasperikova,  2006),  and  the  TD  EMI  array  (Nelson 
and  McDonald,  2001)  have  multiple  transmitters  and/or  receiver  coils.  This  section  describes  in 
detail  the  numerical  procedure  we  use  to  estimate  the  vector  potential  A  and  the  scalar  magnetic 
potential  fusing  multistatic  or  monostatic  EMI  data.  It  is  important  to  notice  that  each  kind  of  data 
requires  a  different  numerical  scheme:  in  the  case  of  multistatic  data  only  one  set  of  equivalent 
responding  magnetic  sources  is  needed  to  predict  an  object’s  EMI  response  at  every  measurement 
point  (Dampney,  1969);  on  the  other  hand,  in  the  monostatic  case  each  data  point  requires  a 
different  set  of  responding  magnetic  sources,  and  this  leads  to  an  underdetermined  ill-posed  least- 
squares  problem.  The  NSMS  model  was  recently  introduced  to  address  these  difficulties.  The 
method  uses  one  set  of  normalized  magnetic  sources  to  predict  monostatic  EMI  data. 

In  this  section  we  present  two  numerical  procedures  to  estimate  the  vector  potential  A  and 
the  scalar  magnetic  potential  y/  starting  from  monostatic  and  multistatic  data.  For  multistatic  data 
we  use  the  Method  of  Auxiliary  Sources  (Shubitidze  et  al ,  2002)  and  for  monostatic  data  we 
employ  the  2D  NSMS  model. 
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Figure  2.2.  Multistatic  EMI  data  survey  scenario  showing  the  intermediate  NSMS 
surface.  The  potentials  A  and  y/  are  obtained  from  sources  on  this 
surface  which  reproduce  the  dipolar  target. 
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2.4.1.  Determining  vector  and  scalar  potentials  from  multistatic  data 

Let  an  EMI  sensor  with  one  or  more  transmitters  excite  a  highly  conducting  and  permeable 
metallic  object  (see  Figure  2.2).  The  object  responds  by  producing  a  secondary  magnetic  field 
which  is  recorded  by  one  or  more  receivers  at  different  locations.  Under  this  scenario,  the  primary 
magnetic  field  induces  only  one  set  of  eddy  currents  or  dipoles  inside  the  object.  Our  goal  is  to 
estimate  the  scattered  magnetic  vector  A  and  scalar  y/  potentials  from  the  measured  multistatic 
magnetic  field  data  and  then  use  these  estimated  potentials,  together  with  the  scattered  magnetic 
field  H,  to  determine  the  buried  object’s  location,  orientation  and  magnetic  polarizability  using 
(2.7)  and  (2.8),  or  (2.11)  and  (2.13). 

First,  choose  a  fictitious  planar  surface  S  underneath  the  set  of  measurement  positions  but 
above  the  target  itself,  and  assume  that  the  scattered  magnetic  field  H"c  due  to  the  target  and  its 
corresponding  vector  and  scalar  potentials  A*c  and  y/T  are  approximated  by  finding  a  set  of 
surface  magnetic  dipoles  placed  on  this  surface  (see  Figure  2.2)  which  produce  a  magnetic  field  that 
most  closely  matches  the  secondary  magnetic  field  produced  by  the  target  itself  at  the  measurement 

locations.  Let  these  approximated  values  be  H^,  A*c ,  and^sc,  where  the  tilde  denotes  that  the 
values  are  approximated  using  the  intermediate  fictitious  surface: 


r  1  f  3R(Rm)  ^  , 

(2.14) 

=-j2*-^G<Rids' 
s  K 

(2.15) 

r(mxR) 

Ar  =  JL-^rJG(R)ds'  . 

(2.16) 

s 


The  auxiliary  (or  fictitious)  surface  S  is  divided  into^,  n=  1,  2,  . . . ,  TV,  subsurfaces,  and  on  it  we 
define  a  local  Cartesian  coordinate  system  using  x,  y,  and  z  (see  Figure  2.2).  According  to 
Huygens’s  principle,  the  scattered  magnetic  field  and  the  vector  and  scalar  potentials  can  be 
reproduced  by  using  a  set  of  surface  magnetic  dipoles  with  only  two  independent  components  along 
x  and  y  .  Thus  the  unknown  surface  dipole  moment  for  layer  S ,  m(s') ,  can  be  written  as 

m (s')  =  Px(s')x  +  Py(s')  y;  (2.17) 

here  Px (s ')  and  Py(s')  are  unknown.  To  determine  P^  (s')  (g  =  {x,y})  for  the  entire  surface  s'  we 
expand  P^ (s')  as 


n= 1 


(2.18) 


where  Pf ,  n  =  1,  2,  3,  ...,  A,  are  expansion  coefficients.  For  computational  simplicity,  in  all 
subsequent  analysis  we  will  assume  the  expansion  functions,  pf(U),  are  a  set  of  pulse  functions 
given  by 
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(2.19) 


FHs1) 


J1’ 

[O,  otherwise. 


This  expansion  of  Fff  (s ?)  results  in  a  piecewise  constant  approximation  to  Px  (, s ?)  and  Py  (s ?)  on 
the  fictitious  surface  S. 

Substituting  (2.18)  into  expression  (2.14)  for  the  magnetic  field  we  obtain 


H 


V)  =  Z  Z  (3R„(R„  ■?&„)- RlL)  ds', 

%=p,z  n= 1  s„  ^ n 


(2.20) 


where  i?w=|r-rw|  is  the  distance  between  the  observation  point  and  the  subsurface  Sn.  The 

potentials  A*c  and  ^sc  can  be  expressed  in  the  same  form  as  (2.20)  after  substituting  (2.18)  into 
(2.16)  and  (2.15)  respectively.  Using  (2.19)  in  (2.20)  enables  us  to  write 


Hr(r)=  X 


±Pf  f  (3R„(R„-i)-i?„2i) 


GiRj 

R: 


ds'. 


(2.21) 


For  convenience,  we  let 


(2.22) 


Then  substituting  (2.22)  in  (2.21)  yields  a  numerical  expression  for  the  magnetic  field  at  r  due  to  a 
magnetic  dipole  layer  m(s') : 


H,c  (r)  =  Pffir,  x, )  +  P2x/(r,  x2 )  +  ...  +  p;/(r,  xw )  + 

^yf(r,y1)  +  P/f(r,y2)  +  ...  +  P^f(r,yN). 

The  physical  interpretation  of  this  equation  is  as  follows.  A  fictitious  surface  located  in  between 
a  target  and  the  measurement  location(s)  consisting  of  a  surface  magnetic  dipole  layer  has  been 
used  to  reproduce  the  magnetic  field  from  an  actual  target.  The  fictitious  surface  has  been  divided 
up  to  N  subsurfaces,  with  the  P;  being  an  unknown  constant  over  each  subsurface.  The  amplitudes 
of  the  Pf  surface  magnetic  sources  can  be  determined  by  matching  the  sum  of  the  magnetic  fields 
produced  from  all  N  subsurfaces  to  the  measured  field  Hdata  (r)  at  r.  Usually,  data  are  recorded  at 
multiple  points  r ,  j  =  1,2,  ...,  J.  After  matching  the  magnetic  field  produced  by  the  fictitious 

surface  layer  to  the  magnetic  field  data  at  J  points,  the  problem  reduces  to  a  linear  system  of 
equations  which  in  compact  matrix  form  can  be  written  as 

[Zy„][^]  =  [Hda,a(ry)]  (2.24) 


where  [Zjn  ]  =  /(r. ,  %n )  is  given  by  (2.22)  and  [Hdata(r)]  is  the  measured  magnetic  field  at  r 
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Using  standard  matrix  algebra  solution  techniques  for  non-square  systems  of  equations,  the 
solution  to  (2.24)  can  be  written  symbolically  as 


[z;„]r[H"-(r;)] 

ixm.] 


(2.25) 


In  order  to  find  a  unique  solution,  J  must  be  larger  than  N.  The  number  of  subsurfaces  N  may  be  as 
few  as  9  resulting  from  a  3  x  3  grid. 

Once  [Pf  ]  is  determined,  all  three  global  values:  estimates  of  the  H  magnetic  field  vector, 

and  the  A  and  ^potentials,  H,c ,  AjC ,  and^sc ,  can  be  computed  readily  at  any  point  in  space  using 
(2. 14)— (2. 16).  After  finding  these  estimates  found  from  monostatic  data,  one  can  find  the  location  R 
and  dipole  moment,  m  (and  thereby  the  orientation,  see  below),  using  (2.7)  and  (2.8)  respectively. 

Note  that  each  subsurface  source  produces  a  separate  Hf ,  A*c ,  and  ^sc  (and  thereby  R  and  m), 
and  in  practice  these  are  averaged  over  N  to  produce  a  final  estimated  R  and  m.  This  concludes  our 
treatment  of  multistatic  data.  In  the  next  subsection  we  estimate  the  location  and  orientation  of  a 
subsurface  target  starting  from  monostatic  measurements. 

2.4.2.  Determining  A*c  and  ^sc  from  monostatic  data  using  NSMS 

Many  currently  available  EMI  sensors,  such  as  the  Geonics  EM-61  and  EM-63  in  the  time 
domain  and  the  Geophex  GEM-3  in  the  frequency  domain,  are  monostatic  and  have  collocated 
transmitters  and  receivers.  Usually  during  UXO  surveying,  the  monostatic  sensor  is  placed  at 
multiple  locations  in  order  to  obtain  more  information  about  a  buried  object:  the  impinging  primary 
magnetic  field  around  the  scatterer  changes  with  location,  and  so  do  the  induced  eddy  currents 
inside  the  object.  Thus  there  is  not  a  unique  set  of  eddy  currents  that  will  correctly  reproduce  a 
given  object’s  monostatic  data  at  multiple  locations  simultaneously.  To  relate  monostatic  data  at 
multiple  locations  with  some  set  of  unique  parameters,  recently  a  new  numerical  technique  called 
the  normalized  surface  magnetic  charge  model  (NSMC)  was  developed  for  a  finite-sized  object 
(Shamatava  et  al.,  2004,  Shubitidze  et  al.,  2007).  Here,  the  NSMC  is  generalized  to  the  NSMS 
model  (by  utilizing  dipole  sources  instead  of  charge  sources)  for  predicting  the  scattered  magnetic 
field  vector  and  the  scalar  potentials  using  only  measured  EMI  magnetic  field  data. 

The  NSMS  approach  assumes  that  the  scattered  magnetic  field  is  produced  by  a  set  of 
magnetic  dipoles  placed  on  a  fictitious  surface  S ,  just  as  in  the  multistatic  case  described  above.  The 
fictitious  surface  is  again  divided  into  small  subsurfaces  and  at  each  subsurface  (and  for  each 
monostatic  measurement)  the  impinging  (primary)  magnetic  field  is  determined.  The  amplitudes  of 
the  responding  NSMS  magnetic  dipoles  at  each  subsurface  are  scaled  by  the  primary  magnetic 
field: 
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Sensor 


Observation  point 


Figure  2.3.  Geometry  of  the  NSMS  model. 


m (V)  =  Px(s')Hpr  x  +  Py(s')Hpr  y. 


(2.26) 


This  is  the  main  difference  between  the  multistatic  section  above  and  the  present  monostatic  NSMS 
formulation.  The  scattered  magnetic  field  in  the  monostatic  case  is  then  determined  to  be 

HT(r)=  Z  £  P:  f  H l  (3R„(R„  ds'.  (2.27) 


£=P>Z  n=l 


After  following  the  same  procedure  as  in  the  previous  section,  the  amplitudes  of  the  NSMC  are 
determined  by  matching  H*c  from  (2.27)  to  the  measured  field,  Hdata(ry),  at  a  selected  set  of 

measurement  points  r hj  =  1,  2,  ...,  J.  Once  the  amplitudes  [Pf]  of  this  source  set  are  found,  the 
scattered  EMI  field  and  its  corresponding  potentials  can  be  determined  at  any  point.  For  more 
details  on  NSMS  refer  to  the  next  section.  Note  that  in  order  to  avoid  ill-conditioning,  the  number 
of  data  points  ./must  not  be  smaller  than  the  number  of  patches  N. 
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2.5.  The  Normalized  Surface  Magnetic  Source  model 


The  NSMS  model  can  be  considered  as  an  extended  version  of  the  multi-dipole  model. 
According  to  the  NSMS  model,  an  object’s  response  to  a  sensor’s  primary  field  is  modeled 
mathematically  by  a  set  of  equivalent  magnetic  dipoles  distributed  over  a  surface  surrounding  the 
object.  The  amplitudes  of  the  equivalent  responding  magnetic  dipoles,  which  are  oriented  normally 
to  the  surrounding  surface,  are  normalized  by  the  component  of  the  primary  magnetic  field  normal 
to  the  surface.  The  amplitude  Q  of  the  dipole  sources  is  a  property  of  the  object,  and  its  integral  Q 
over  the  surface  constitutes  a  global  magnetic  “capacitance”  of  sorts  that  relates  the  induced  double 
layer  charge  density  to  the  normal  component  of  the  primary  magnetic  field,  //f .  Thus  the 
scattered  magnetic  field  of  the  NSMS  is  identical  to  that  produced  by  a  layer  of  magnetic  dipoles 
distributed  over  the  surface.  It  can  be  considered  that  such  a  surface  distribution  is  generated  by 
spreading  positive  charge  density  over  the  positive  side  (with  respect  to  the  surface  normal)  of  a 
smooth  surface  and  an  identical  distribution  of  opposite  sign  on  its  negative  side.  The  result  is  a 
double  layer  of  magnetic  charge  separated  by  an  infinitesimal  distance.  According  to  the  Huygens 
equivalence  principle,  an  object’s  EMI  response  can  be  reproduced  with  a  set  of  equivalent  single 
or  double  layer  magnetic  charges.  The  double  layer  charge  density  introduces  the  proper 
discontinuities  in  the  tangential  components  of  magnetic  flux  density  vector  B,  but  does  not  affect 
the  transition  of  the  normal  component  of  B.  The  single  layer  charge  density,  on  the  other  hand,  in 
no  way  affects  the  transition  of  the  tangential  components  but  accounts  for  the  proper  discontinuity 
in  the  normal  component  of  the  vector  B.  Since  the  normal  component  of  the  magnetic  flux  density 
vector  is  always  continuous  on  the  boundary  between  two  media  (there  are  no  free  magnetic 
charges  in  nature),  the  secondary  magnetic  field  outside  the  object  can  be  accounted  for  by  a  layer 
of  magnetic  dipoles  with  density  Q  distributed  on  an  auxiliary/surrounding  surface. 

In  order  to  illustrate  the  similarity  between  the  multi  dipole  and  NSMS  approaches,  let  us 
briefly  describe  the  multi-dipole  model  approximation  (Shubitidze  et  al .,  2008).  In  the  multi  dipole 
approximation  the  EMI  response  of  the  whole  object  is  represented  by  a  set  of  magnetic  dipoles 
distributed  on  a  line  along  the  object’s  elongation  axis  (red  line  on  Figure  2.3).  Mathematically,  the 
secondary  magnetic  field  due  to  a  set  of  dipoles  of  moment  m(f ?)  along  the  line  L  an  be  written  as 


Hse  =  [ - - — -m(f')-(3RR- 1)  dT,  (2.28) 

{ 

where  R  is  the  unit  vector  along  R  =  r  —  r  ,  r  is  the  position  of  the  dipole  on  line  L,  r  is  the 
observation  point  (see  Figure  2.1),  I  is  the  identity  dyad,  and  m(f ')  is  the  induced  magnetic  dipole 
along  L.  The  induced  magnetic  dipole  is  related  to  the  primary  field  via  m(f')  =  M(f')-H 
where  the  polarizability  tensor  \1  ( i ')  is  assumed  to  be  varying  along  the  line  and  to  have  the  form 


M(f) 


mjn  o  0 
0  mjn  0 
0  0  mzAJ:') 


(2.29) 
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where  mxx{V) ,  myv (£ ') ,  and  mzz(£')  are  the  moments  along  the  principal  axes  of  M(f ?)  at  point  V 

on  L.  The  multi  dipole  model  is  a  fast  algorithm,  and  it  takes  into  account  an  object’s  geometric  and 
electromagnetic  (heterogeneity)  variations  along  L.  However,  actual  objects  have  finite  size,  and  the 
primary  magnetic  field  can  vary  not  only  along  the  axis  of  symmetry  but  along  the  circumference 
axis  as  well.  To  take  into  account  geometric  variations  along  the  circumference  axis  let  us  extend 
the  multi  dipole  model  by  distributing  the  dipoles  on  a  closed  surface.  Thus,  the  total  scattered 
magnetic  field  produced  by  the  surface  magnetic  dipoles  (i.e.,  the  NSMS  system)  can  be  written  as 

h  c = f  ■  <3M  -  b ds'  ’  (2-3°) 

s  4/riU»R 

where  the  induced  magnetic  dipole  m(r  ,)  on  the  surface  at  point  r,  is  related  to  the  primary  field 
via  m(r  ,)  =  M(rv)-H;?r  (r ,).  In  the  NSMS  approach  the  polarizability  tensor  M(rs,)  in  a  local 
curvilinear  coordinate  system  (u(r  ,) ,  v(r  ,) ,  n(r5,) )  is  assumed  to  have  only  one  nonzero  principal 
axis  along  the  normal  direction  n(rs.) ,  and  since 

H"  (O  =  K (r,)u(r s,)  +  Hr  (r,)v(r ,)  +  H?  {r,)n(r,) 


we  have 


m(rs,)  =  Q.(rs,)H„r  (rs,)n(rs,) ,  (2.31) 

where  Q(r5,)  is  the  amplitude  of  the  magnetic  dipole  at  point  r  ,  on  the  spheroidal  surface.  Overall, 
in  the  NSMS  approach  the  scattered  magnetic  field  is  approximated  as  a  sum  of  magnetic  fields 
radiated  by  elementary  magnetic  sources  (magnetic  charges  or  dipoles)  distributed  on  an  auxiliary 
closed  surface.  The  amplitudes  of  these  responding  magnetic  charges  at  each  point  on  the  surface 
are  normalized  by  the  normal  component  of  the  impinging  primary  magnetic  field,  .  At  the  end 
the  scattered  field  is  expressed  by  the  proportionality  constant  Q(sf) : 

HsVr)=j— (rs.)  (n(rs,)-(3M-T))  ds'  .  (2.32) 

s  ^0K 

The  quantity  essentially  a  normalized  charge  density,  corresponding  to  a  primary  field 

H„r  =  1  A/m.  For  discrimination  purposes  we  propose  using  a  single  variable  equal  to  the  total 

normalized  magnetic  charge.  This  is  a  global  measure  of  Q  for  the  entire  object  and  is  thus  less 
subject  to  numerical  fluctuation  than  the  individual  Cl(t,s')  values: 

Q(t)  =  jn(t,s')ds'  (2.33) 

S 

Our  studies  have  showed  that  the  total  NSMS  Q  is  invariant  for  a  given  object,  in  the  sense 
that  different  computational  constructs  ( e.g .,  surrounding  surfaces)  and  different  primary  fields 
produce  the  same  value.  Once  the  amplitudes  of  the  NSMS  are  determined  the  scattered  EM  fields 
can  be  calculated  extremely  fast  and  accurately  for  any  known  object  and  any  type  of  EMI  sensor. 
Thus  the  NSMS  model  can  be  considered  a  complete,  unified  fast  forward  dipole  model  of  EMI 
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scattering  for  use  in  discrimination  processing  that  relies  on  comparing  received  fields  to  those 
produced  by  a  library  of  known  objects.  In  such  “pattern  matching”  discrimination  processing,  an 
optimization  or  search  algorithm  is  used  to  determine  what  configuration  of  known  UXO  types  and 
sensor  best  reproduces  the  received  signal,  for  catalogued  Q(t,sf)  distributions  associated  with 
those  UXO. 


2.6.  A  numerical  procedure  for  estimating  orientation 

The  Method  of  Auxiliary  Sources  (MAS)  (Shubitidze  et  al.  2002),  approximates  an  object’s 
response  to  a  given  excitation  as  the  summation  of  magnetic  fields  produced  by  elementary 
magnetic  dipoles  or  charges  placed  inside  the  object.  Using  the  superposition  principle,  this  set  of 
dipoles  can  be  further  approximated  as  one  independent  aggregate  dipole  which  can  then  be 
compared  to  the  simple  dipole  model  commonly  in  use.  According  to  the  simple  dipole  model 
(Pasion  and  Oldenburg,  2001),  the  secondary  magnetic  field  due  to  the  dipole  m  is 

H  =  — !-_m.(3RR-T),  (2.34) 

4  KjU0R 

where  R  is  the  unit  vector  along  U  r  r  .  r  is  the  dipole’s  position,  r  is  the  observation  point, 

and  I  is  the  identity  dyad  (see  Figure  2.1).  The  dipole  moment  m  induced  by  the  primary  magnetic 
field  Hpr  is  given  by 


m  =  M  •  Hpr  (r,r^) , 


(2.35) 


where  M  is  the  target’s  magnetic  polarizability  tensor.  This  tensor  depends  on  the  scatterer’s 
shape,  size,  and  material  properties.  In  a  coordinate  system  not  aligned  with  the  scatterer’s  principal 
axes  for  different  primary  magnetic  fields  Hpr(ry,r^)  ,y  =  1,2,  . ..,  J,  equation  (2.35)  can  be  written 
in  matrix  form  as 


where 


[m]  =  M.[H^], 


(2.36) 


m'x  ,... 

■«T(r.,r)  .... 

.  f 

m\  ,.. 

,  mi 

and 

[H/,']= 

«f(r,  ,r)  ... 

•>  7/f(r„r) 

m\  ,... 

,  mJz\ 

1 

■'i 

A 

H?{r„  r) 

(2.37) 


Equation  (2.36)  can  be  solved  for  the  magnetic  polarizability  M : 


[H'”  ]  f 


(2.38) 
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Once  M  is  determined,  the  principal  target  and  the  global  coordinate  systems  can  be 
related  via  the  Euler  rotation  tensor,  A(i//,0,</>) ,  which  is  an  orthogonal  matrix:  A-1  =  Ar  .  Using 
the  eigendecomposition  theorem  one  can  write  M  as 


M  =  A  J3  A\ 


where 


P  = 


0  fi„ 
0  0 


0 

0 


P=\ 


(2.39) 


(2.40) 


For  a  body  of  revolution  (BOR)  (which  includes  most  UXO)  fixx  =  /3yy,  and  because  of  BOR 
symmetry  the  third  Euler  angle  y/  is  zero,  and  thus 


A  = 


cos  # cos (j)  cos  # sin (j)  -sin# 
-  sin  (j>  cos  (j)  0 

sin#cos^  sin#sin0  cos# 


(2.41) 


where  #and  (f)  are  the  angles  between  the  local  and  global  axes.  Equation  (2.39)  can  be  rewritten  as 


at  ma  =  p . 


(2.42) 


Since  M  is  determined  from  (2.38)  and  the  A  are  given  explicitly,  from  (2.42)  we  can  determine  6 
and  Rvalues  that  give  minimum/zero  values  of  the  off-diagonal  elements  of  /?  . 


2.7.  Matrix  diagonalization 

To  determine  the  matrices  A  and  /?  we  employ  a  simultaneous  matrix  diagonalization 
technique  that  uses  the  Jacobi  angles  to  build  a  unitary  orthogonal  matrix  A  .  Namely,  consider  a  set 
|  M  ={Mt  1 1  =  1,  Af }  of  A,  3  x  3  matrices.  When  the  matrices  in  |M  are  normal  commuting 
matrices  ,  their  off-diagonal  terms  can  be  set  to  zero  by  a  unitary  transform,  thus  simultaneously 
diagonalizing  the  set  |  M  . 


min 


(of  |  M)  =  £ 

\<it *  jt<  3 


M 


it, ft 


(2.43) 


M ;!  denotes  the  (it,  /f)-th  entry  of  matrix  |  M  .  The  simultaneous  diagonalization  is  obtained  by 

minimizing  the  composite  objective  function  ^  |  ,  l\T.  l7  |2  by  a  unitary  matrix  A  . 
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2.8.  Results 


2.8.1.  Estimating  the  location  of  a  buried  object  using  synthetic  EMI  data 

In  this  section  we  present  results  of  numerical  experiments  to  demonstrate  the  applicability 
of  the  physics-based  method  for  determining  a  buried  object’s  location  without  solving  a  time- 
consuming  minimization  problem.  First,  the  proposed  algorithm  is  tested  for  the  case  of  a  single 
dipole  source,  with  dipole  moment  m,  placed  in  free  space  (k  =  0)  at  different  depths  below  an 
arbitrary  measurement  plane  with  no  primary  field.  In  these  tests,  the  measurement  plane  was 
divided  into  grid  points  (40  x  40),  and  at  each  point  on  the  grid  the  magnetic  field,  scalar  and  vector 
potentials  were  calculated  from  the  single  dipole  source  on  the  fictitious  surface  using  equations 
(2.1)-(2.3)  for  the  magnetic  dipole  m  =  (1,1,2).  Then,  using  these  values  of  H,  A,  and  y/,  the 
dipole’s  location  was  determined  using  (2.7)  at  each  grid  point.  The  synthetic  test  is  done  for  four 
cases:  h=  1,  10,  30,  and  40  cm  below  the  measurement  surface.  The  inverted  location  of  the 
magnetic  dipole  is  depicted  in  Figure  2.4. 

The  results  show  that  the  proposed  algorithm  inverts  the  dipole’s  location  from  each  point  very 
accurately,  i.e.  for  determining  the  dipole’s  location,  only  H,  A,  and  y/  at  a  single  point  (each  grid 
point  separately)  are  needed.  In  addition,  results  indicate  that  the  approach  is  applicable  in  the  near 
(h=\  cm)  as  well  as  far  field  (h  =  41  cm). 

In  order  to  test  the  effects  that  the  presence  of  multiple  dipoles  would  have  on  the  algorithm 
(which  assumes  only  one  dipole  source),  similar  tests  were  done  for  three  separate 


source  dipoles,  each  with  m  =  (1,1,2),  placed  at  [x,  y,  z\  =  [-5,  0,  14]  cm,  [x,  y,  z]  =[-5,  0,  19]  cm 


Figure  2.4.  Inverted  coordinates  (a)  x0,  (b)  y0  and  (c)  z0  for  a  dipole.  Inverted  locations  for 
these  synthetic  tests  are  essentially  exact  and  demonstrate  the  viability  of  this 
method. 


Figure  2.5.  Three  source  dipoles  at  different  depths:  (a)  Problem  geometry.  Inverted  coordinates 
of  inverted  aggregate  dipole  (b)  x09  and  (c)  z0. 

2.8.2  Predicting  a  dipole’s  orientation  from  monostatic  synthetic  data 

Here  we  test  the  new  approach  for  estimating  a  dipole’s  orientation  from  monostatic  field 
data.  The  dipole  with  magnetic  polarizability  /?  tensor  diag (/?  )  =  (1,  2,  3)  is  placed  30  cm  bellow 
the  measurement  surface.  The  dipole  is  illuminated  which  the  virtual  EM-63  sensor,  and  mono¬ 
static  sensor  data  are  collected  over  a  60  cm  x  60  cm  measurement  surface  divided  into  25  x  25  grid 
points.  The  mono-static  means  that  the  transmitter  and  receiver  are  always  co-located  i.e.  when  it  is 
located  at  some  particular  point  on  measurement  surface  the  sensor  records  only  the  scattered  field 
at  that  point  produced  by  the  primary  field  transmitted  from  that  point.  Each  data  point  is  from  a 
different  scattered  field,  as  each  result  from  a  different  excitation  of  the  dipole.  At  each  point  on  the 
grid  the  magnetic  field,  scalar  and  vector  potentials  were  calculated  from  the  magnetic  dipole  with 

polarizability  (3  tensor  diag(/?  )=(1,  2,3)  on  the  measurement  surface  using  (2. 1)— (2.3)  and  the 

simple  dipole  model  approximation  (2.35).  First,  using  these  values  of  H,  A,  and  y/9  the  dipole’s 
location  was  estimated  using  (2.7),  and  then  using  the  predicted  values  of  location,  the  dipole’s 

magnetic  polarizability  tensor  M  was  determined  from  (2.8)  and  (2.38).  Finally,  the  dipole’s 
orientation  (y/9  6 \  (/))  that  gives  minimum/zero  values  of  the  /?  matrix’s  off-diagonal  elements  and 

diagonal  elements  of  [3  were  estimated.  The  numerical  tests  are  done  for  five  cases  and 
summarized  in  Table  2.1.  For  these  numerical  studies,  a  comparison  between  estimated  and  true 

values  for  the  orientation  shows  an  excellent  match.  Note  that  the  estimated  values  of  (3  diagonal 
elements  exactly  matched  with  the  actual  diag(/?  )  =  (1,  2,  3)  values  as  well. 


2.8.3  Noise  effect:  predicting  a  dipole’s  orientation  from  monostatic  noisy  synthetic  data 

In  this  section  we  studied  the  robustness  of  the  reduced  HAP  approach  with  respect  to  noise 
for  a  dipole  source.  We  consider  the  same  dipole  with  magnetic  polarizability  /?  tensor 

diag(/? )  =  (1,  2,  3)  oriented  with  Euler  angles  (0=  45°,  </>=  30°,  y/=  60°)  and  placed  35  cm  below 
the  measurement  surface.  The  dipole  is  illuminated  with  the  virtual  EM-63  sensor,  and  monostatic 
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sensor  data  are  collected  over  a  40  cm  x  40  cm  measurement  surface  divided  into  21x21  grid 
points.  The  distribution  of  the  scattered  magnetic  field  and  the  scalar  potential  without  (left)  and 
with  (right)  noise  are  depicted  on  Figure  2.6.  As  before  at  each  point  on  the  grid  the  magnetic  field 
and  scalar  potential  were  calculated  from  the  magnetic  dipole,  on  the  measurement  surface  using 
(2. 1)— (2.3)  via  the  basic  simple  dipole  model  approximation  (see  distribution  on  the  left  column  of 
Figure  2.6).  Then  at  each  point  on  the  measurement  surface  a  random  Gaussian  random  noise  were 
to  both  field  and  scalar  potentials.  The  noise  generated  this  way  is  better  suited  to  simulating  sensor 
noise  than  to  modeling  positional  uncertainty.  The  noisy  data  distribution  are  plotted  on  Figure  2.6 
left.  Using  these  values  of  H  and  ^(Figure  2.6,  left),  the  dipole’s  location  was  estimated  by  using 
(2.11)  and  then,  using  the  predicted  values  of  location,  the  dipole’s  magnetic  polarizability  tensor 

M  was  determined  from  (2.13)  and  (2.38).  Finally,  using  the  joint  orthogonal  diagonalization 

technique  (2.43)  the  M  matrix’s  diagonal  elements  and  orientations  (<9,  y/)  were  estimated.  The 

numerical  test  is  summarized  in  Table  2.2.  For  these  numerical  studies,  the  comparison  between 
estimated  and  true  values  for  the  orientation  show  very  close  matches. 


Table  2.1.  Estimated  and  actual  orientation  of  a  dipole  using  a  synthetic  data  set. 


Estimated 

Actual 

Case 

9  (°) 

4>(°) 

'p  (°) 

e  (°> 

*(°) 

'p  (°> 

1 

45 

3 

1.8 

45 

3 

1.8 

2 

45 

30 

18 

45 

30 

18 

3 

60 

45 

90 

60 

45 

90 

4 

30 

20 

60 

30 

20 

60 

5 

36 

45 

45 

36 

45 

45 

Table  2.2.  Estimated  and  actual  orientation  for  the  dipole  using  noisy  data. 


Estimated 

Actual 

X0(cm)  Y0(cm)  Z0(cm)  0(0)  <t>(°)  Y  (°)  Pi  P2  P3 

X0  (cm)  Y0(cm)  Z0  (cm)  0(0)  <j,(°)  v  (0)  Pl  p2  p3 

2.2  1.5  -35.36  44  34  64  1.03  1.98  2.96 

0  0  -36  45  35  60  1  2  3 
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Figure  2.6.  Scattered  magnetic  field  and  scalar  magnetic  potential  for  a  dipole  source.  Left:  no 
noise  added;  right:  noise  added. 
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Figure  2.7.  Problem  geometry  for  a 
sphere  under  the  meas¬ 
urement  plane. 
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Figure  2.8.  Estimated  zc  for  the 
sphere  in  Figure  2.7. 
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Figure  2.9.  Estimated  yj x 2C  +  y2c  for 
the  sphere  in  Figure  2.7. 


Figure  2.10.  Estimated  zc  for  the 
sphere  in  Figure  2.7 
vs.  observation  point. 
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Figure  2.11.  Geometry  of  the  problem 
for  a  homogeneous  object 
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Figure  2.13.  Estimated  ^x2  +  y2  for  the 
bullet  of  Figure  2.11. 
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Figure  2.12.  Estimated  zc  for  the  bullet 
of  Figure  2.11. 
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Figure  2.14.  Estimated  zc  for  the  bullet 
in  Figure  2.1 1  vs. 
observation  point. 
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2.8.4.  Studies  for  a  highly  conducting  and  permeable  sphere 

Next,  the  algorithm  was  applied  to  the  synthetic  EMI  response  from  a  highly  conducting 
and  permeable  metallic  sphere,  with  radius  a  =  5  cm,  permeability  jur=  100,  and  conductivity 
<j=  4  x  106  S/m.  The  sphere  is  placed  at  [xc,yc,  zc ]  =  [0,0-25]  cm  (Figure  2.7)  and  illuminated  with 
a  time-varying  magnetic  field  generated  by  a  virtual  Geophex  GEM-3  sensor  (Won  et  al. ,  1996). 
The  sensor’s  primary  magnetic  field  is  modeled  as  the  field  radiated  by  two  concentric  coils.  The 
current  amplitudes  in  the  coils  are  chosen  in  such  a  way  that  at  the  center  of  the  concentric  coils, 
where  the  receiving  coil  resides,  the  primary  magnetic  field  is  zero.  The  measurement  plane  is  again 
divided  into  grid  points,  and  the  GEM-3  sensor  is  swept  over  the  grid  points  in  25-cm  increments 
resulting  in  a  5  x  5  grid  of  measurement  locations.  At  each  measurement  position  a  separate  EMI 
problem  is  solved  using  the  MAS  EMI-BOR  code  and  the  scattered  magnetic  field  and  scalar  and 
vector  potentials  (Hisc,  Aisc,  and  y/ isc)  are  calculated  at  the  sensor’s  center,  i.e.  a  synthetic  mono¬ 
static  data  set  is  created.  Using  this  Hisc,  Aisc,  and  y/ isc,  the  sphere’s  center  [xc,yc,  zc\  is  determined 

using  (2.7)  at  each  grid  point.  Figures  2.8  and  2.9  show  the  predicted  zc  and  ^x2c  +  y2c  over  the  x-y 

measured  plane  using  the  analytic  inversion  algorithm  presented  in  this  paper.  The  results  show  that 
the  algorithm  predicts  the  sphere’s  center  location  with  reasonable  accuracy  at  each  measurement 
point.  Figure  2.10  shows  a  comparison  between  the  predicted  and  actual  zc  for  each  observation 
point.  The  results  demonstrate  that  knowing  Hisc,  Aisc,  and  y/\c  at  any  given  point  is  enough  to 
estimate  the  buried  sphere’s  location. 


2.8.5.  Studies  for  non-spherical  and  heterogeneous  objects 

The  sphere  is  the  simplest  3D  geometrical  figure,  and  its  EMI  response  can  be  represented 
very  well  with  a  magnetic  dipole.  However,  in  order  to  demonstrate  the  applicability  of  this  new- 
physics  based  algorithm,  in  this  section  we  present  additional  numerical  studies  for  homogeneous 
(Figures  2.11-2.14)  and  heterogeneous  (Figures  2.15-2.18)  objects.  The  homogeneous  object’s 
geometrical  and  electromagnetic  parameters  are:  nose  radius  a  =  2.5  cm,  overall  length  18  cm, 
relative  permeability  jur  =  100,  and  conductivity  cr  =  4xl06  S/m.  The  heterogeneous  object’s 
parameters  are:  overall  length  28  cm,  nose  radius  2.5  cm,  relative  permeability  /jr=  100,  and 
conductivity  o  =  4xl06  S/m  for  the  bottom  cylindrical  part,  while  the  nose  has  /jr  =  30  and 
cr  =  5xl07  S/m. 

The  objects  are  placed  at  \xc,  yC9  zc]  =  [0,  0,  -25]  cm  (see  Figures  2.11  and  2.15)  and 
oriented  at  45°  with  respect  to  the  vertical.  As  in  the  case  of  the  sphere,  they  are  illuminated  with  a 
time-varying  magnetic  field  generated  by  a  synthetic  double-loop  model  of  the  Geophex  GEM-3 
sensor.  The  simulated  sensor  is  placed  at  different  points  on  the  measurement  grid  at  25-cm 
increments.  The  full  EMI  problem  is  solved  for  each  sensor  position,  and  H*c ,  A*c ,  and  y/T  are 
calculated.  The  objects’  center  location  is  estimated  using  (2.7)  and  is  depicted  in  Figures  2.12, 
2.13,  2.16,  and  2.17.  These  results  again  demonstrate  that  the  algorithm  can  predict  a  buried 
object’s  center  position  with  reasonable  accuracy  for  both  homogeneous  and  heterogeneous  objects. 
Note  that  these  results  are  presented  only  for  a  single  frequency;  similar  trends  have  been  observed 
for  all  other  frequencies.  Figures  2.13  and  2.17  show  the  buried  object’s  center  in  the  x-y  plane, 
while  Figures  2.14  and  2.18  depict  the  estimated  zc  from  each  of  the  observation  points. 
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Figure  2.15.  Problem  geometry  for  a 
heterogeneous  object. 
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Figure  2.17.  Estimated  y  x2c  +  y2c  for  the 

heterogeneous  object  of 
Figure  2.15. 
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Figure  2.16.  Estimated  zc  for  the 

heterogeneous  object  of 
Figure  2.15. 


Figure  2.18.  Estimated  zc  for  the  bullet  of 
Figure  2.15  vs.  observation 
point. 
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Figure  2.19.  Scattered  magnetic  field  from  an  81 -mm  UXO  modeled  via  the  SEA  using  synthetic 
data  a)  without  and  b)  with  noise,  c)  Inverted  location  from  noise-free  and  noisy  data. 
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2.8.6.  Noise-effect  analysis  for  heterogeneous  objects 

To  further  illustrate  the  applicability  of  the  HAP  method  in  real  field  conditions,  we 
investigate  the  performance  of  the  technique  for  a  case  in  which  signals  are  contaminated  with 
additive  sensor  noise.  We  generate  synthetic  scattered-field  data  using  the  (spheroidal)  standardized 
excitation  approach  (SEA).  The  analysis  is  done  for  an  81 -mm  UXO  for  which  the  amplitudes  of  a 
reduced  set  of  responding  sources  were  determined  and  stored  in  a  library.  Using  these  sources  we 
generated  the  object’s  EMI  response  on  a  11x11  measurement  grid  to  the  primary  field  of  the 
(monostatic)  EM-63  sensor.  These  responses  are  depicted  on  Figure  2.19,  upper  left;  the  target  was 
placed  at  a  depth  of  35  cm  and  oriented  horizontally.  To  replicate  the  conditions  of  a  EM-63  survey 
we  calculated  only  one  component  (7/z)  of  the  scattered  magnetic  field.  The  2D  NSMS  layer  was 
placed  under  the  measurement  grid  and  the  amplitudes  of  the  NSMS  were  determined,  after  which 
the  magnetic  field  and  scalar  magnetic  potential  were  determined  for  26  time  channels.  The  target 
position  was  estimated  using  the  HAP  technique.  The  inverted  results  are  depicted  on  Figure  2.19. 
The  results  show  that  the  HAP  method  reconstructs  the  position  very  well.  We  then  added  random 
Gaussian  noise  to  the  field  reading  at  every  position;  the  results  appear  in  Figure  2.19  upper  right. 
The  object’s  position  was  inverted  from  the  noisy  data  by  applying  the  same  numerical  procedure. 
The  inverted  results  are  depicted  on  Figure  2.19  with  red  circles.  We  see  that  even  with  noisy  data 
the  method  does  a  very  good  job  of  predicting  the  object’s  position,  in  particular  its  depth. 


2.8.7.  Estimating  object  location  using  measured  monostatic  and  multistatic  EMI  data 


2.8. 7.  a.  Frequency-domain  EMI  data 

In  this  section  we  apply  our  approach  to  monostatic  GEM-3  data.  We  first  compare  the 
predicted  EMI  magnetic  field  computed  via  the  NSMS  fictitious  source  layer  to  measured  data.  In 
these  tests,  the  NSMS  model  is  used  to  predict  the  scattered  magnetic  field  components  Hx  and  Hy 
using  only  monostatic  measured  magnetic  field  data.  The  tests  are  done  for  a  steel  sphere  of  radius 
a  =  3.75  cm. 

The  data  were  collected  at  two  elevations:  h\  =  2\  cm  and  h2  =  26  cm  using  the  new 
Geophex  GEM-3D+  sensor  (O’Neill  et  al. ,  2004).  Here  h\  and  h2  are  measured  from  the  sphere’s 
center  to  the  measurement  surface.  The  GEM-3D+  sensor  consists  of  two  transmitter  loops  and 
three  orthogonal  receivers  (Figure  2.20).  The  currents  in  the  transmitter  loops  circulate  in  opposite 
directions  and  are  scaled  so  that  their  respective  primary  fields  cancel  at  their  common  center, 
where  the  receiving  coils  are  located.  In  the  numerical  model,  the  transmitter  loops  are  idealized  as 
infinitely  thin  loop  sources  of  radii  a\  and  a2 ,  with  currents  I2  =  -I\  a2la\.  At  each  elevation,  the  EMI 
field  is  measured  on  an  81 -point  grid  at  5 -cm  increments. 
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Figure  2.20.  The  GEM-3D+  sensor  head  Figure  2.21.  Comparisons  between  actual 

showing  the  transverse  and  modeled  H_  data, 

receiver  coils. 

The  scattered  magnetic  field  is  modeled  as  a  superposition  of  the  magnetic  field  produced 
by  a  set  of  normalized  magnetic  dipoles  placed  on  a  fictitious  surface  just  below  the  actual 
measurement  grid.  The  fictitious  surface  is  divided  into  small  patches  (cf  Figure  2.2).  At  each  patch 
the  primary  magnetic  field  is  determined.  The  amplitudes  of  the  responding  magnetic  dipoles  on 
each  subsurface  are  scaled  by  the  primary  magnetic  field  using  (2.26).  The  amplitudes  of  the 
responding  surface  layer  NSMS  are  then  determined  by  matching  the  z-component  of  the  modeled 
magnetic  field  to  its  measured  counterpart  at  the  grid  points  with  a  sensor  elevation  h\.  Once  the 

amplitudes  are  found,  the  secondary  field  quantities  H*c ,  A*c ,  and  ^sc  can  be  found  anywhere 
above  the  fictitious  surface.  Note  that  the  auxiliary  NSMS  layer  is  placed  between  the  object  and 
the  measurement  plane.  In  this  arrangement,  the  NSMS  surface  sources  represents  the  magnetic 
fields  radiated  by  all  metallic  objects  distributed  beneath  the  fictitious  surface. 

The  NSMS  sources  were  distributed  on  a  5  x  5  grid  at  an  elevation  h  =  10.5  cm.  The 
spacing  between  sources  in  this  case  was  10  cm.  The  NSMS  modeled  magnetic  field  is  matched  to 
the  measured  data  at  the  elevation  h\  =  21  cm.  The  results  appear  in  Figure  2.21.  The  surface- 
distributed  NSMS  sources  do  a  good  job  of  predicting  the  EM  signal  at  h\.  After  determining  the 
amplitudes  of  the  NSMS  we  can  compute  the  x-  and  y-components  of  the  scattered  monostatic 
magnetic  field  at  9  x  9  grid  points  at  the  elevation  h2=26  cm,  where  the  GEM-3D+  data  were 
collected.  The  results  appear  in  Figure  2.22.  The  NSMS  dipoles  can  be  seen  to  predict  Hx  and  Hy  of 
the  scattered  magnetic  field  at  elevation  h2  using  only  the  Hz  of  the  scattered  magnetic  field  data  at 

h. 
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Figure  2.22.  Comparisons  between  actual  and  predicted  data  for  a  sphere  subject  to  GEM-3D+ 
excitation:  (a)  x-component,  (b)  ^-component. 


Similarly,  using  these  amplitudes  of  the  NSMS  layer  we  estimate  the  vector  and  scalar 
potentials  AjC  and  ^sc  at  the  same  9x9  grid  at  elevation  h2  =  26  cm  in  order  to  determine  the 

coordinates  of  the  center  of  the  sphere  using  (2.7).  The  predicted  +  y]  and  zc  are  depicted  in 

Figure  2.23.  The  results  illustrate  that  the  algorithm  can  predict  the  object’s  location  from  actual  Hz 
data. 


To  demonstrate  the  superior  performance  of  the  HAP  technique  combined  with  the  NSMS 
approach,  we  apply  the  algorithm  to  frequency-domain  GEM-3 D+  blind-test  data.  These  were 
collected  for  eight  targets,  of  which  three  were  clutter  and  five  were  UXO:  a  BLU-26  bomblet,  a 
Rockeye,  an  81 -mm  projectile,  a  60-mm  mortar,  and  a  57-mm  bullet.  For  all  these  items  detailed 
data  were  collected  under  controlled  conditions,  which  were  used  to  determine  the  total  NSMS  Q(f) 
for  each  target.  We  then  determined  the  locations  and  orientations  of  the  objects  using  the  HAP 
algorithm.  The  true  and  estimated  locations  and  orientations  of  the  targets  are  given  in  Table  2.3. 
The  results  show  that  inverted  and  true  values  of  position  and  orientation  are  in  general  very  close. 
After  that  we  calculated  the  total  NSMS  for  each  of  the  eight  objects  and  compared  them  to  the 
total  NSMS  values  stored  in  a  library  of  six  UXO  items.  In  all  cases  the  inverted  NSMS  closely 
matches  its  corresponding  UXO  items  (see  Figure  2.24),  which  means  that  the  algorithm  was  able 
to  identify  all  items  with  100%  accuracy. 
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Figure  2.23.  Estimated 


(left)  and  zc  (right)  for  the  steel  sphere.  The  true  zc  =  26  cm. 


Table  2.3.  Estimated  and  actual  positions  and  orientations  for  the  GEM-3D+  blind-test  data  runs. 


Blind 

Estimated  (Actual) 

test  # 

0(°) 

x0  (cm) 

yo  (cm) 

z0  (cm) 

1 

93.60(45) 

90.0  ( 90) 

3.0  (  0) 

0.01(  1.5) 

-29.68(-28) 

2 

86.17  (90) 

20.9  (135) 

-5.5  (  -5) 

5.03(  5) 

-34.46(-31) 

3 

105.36  (  0) 

76.4  (arb) 

-15.6  (-20) 

-15.8  (-15) 

-43.46(-42) 

4 

91.08  (arb) 

2.25(arb) 

18.2(20) 

-0.5  (  0) 

-25. 1 9(-22) 

5 

81.05  (135) 

10.11(135) 

20.2  ( 15) 

-0.7  (  5) 

-26.22(-26) 

6 

86.34  (180) 

67.05(arb) 

5.3  ( 10) 

-12.64(-20) 

-20.25(-18) 

7 

65.56  ( 90) 

16.63(45) 

12.9  ( 18) 

4.13(  12) 

-25.81  (-20) 

8 

66.06  ( 90) 

14.01(90) 

3.35  (  5) 

15.45(  10) 

-43.10(-40) 
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Figure  2.24.  Comparison  between  GEM-3D+  blind  tests  inversions  and  library  NSMS. 


2.8.7.b.  Time-domain  EMI-63  data 


The  algorithm  was  also  applied  to  data  from  the  time  domain  EM-63  sensor  made  by 
Geonics  LTD  (McNeill  and  Bosnar,  1996).  The  EM-63  sensor  consists  of  a  1  m  x  1  m  square 
transmitter  loop  and  two  receiver  loops:  (1)  the  main  receiver  loop  is0.5mx0.5min  size  with  a 
center  that  coincides  with  the  transmitter  coil’s  center,  and  (2)  the  same  size  receiver  loop  that  is 
placed  60  cm  above  the  main  receiver  coil.  These  receivers  measure  the  complete  transient  response 
from  180  ps  to  25  ms.  The  time-domain  measurements  were  conducted  at  the  UXO  test-stand  at  the 
U.S.  Army  Corps  of  Engineers,  Engineer  Research  and  Development  Center  (ERDC)  Laboratory, 
Vicksburg,  MS.  The  measurement  platform,  which  is  made  with  nonmetallic  fiberglass  materials, 
has  a  usable  measurement  area  of  3  x  4  m  in  extent.  The  sensor  is  mounted  on  a  robotic  arm  that 
can  be  moved  and  controlled  around  the  measurement  area  using  software  running  on  a  remote 
personal  computer.  Automated 
remote  controls  are  used  to  position 
an  ordnance  item  at  an  accurate 
depth  (within  1  cm)  below  the 
measurement  area.  The  sensor  can 
be  positioned  with  an  accuracy  of 
approximately  1  mm.  All  data 
presented  here  were  collected  by 
Sky  Research  personnel. 

The  time  domain  data  were 
collected  at  a  h  =  60  cm  elevation  at 
441  points  for  a  horizontally 
oriented  81 -mm  UXO.  Again,  the 
scattered  magnetic  field  is  modeled 
as  a  superposition  of  the  magnetic 
field  produced  by  a  set  of  surface 
magnetic  dipoles  placed  on  a 
fictitious  surface  just  20  cm  below 
the  real  surface.  After  determining 
the  amplitudes  of  the  NSMS  on  this 
surface,  the  magnetic  field  and  the 
potentials  are  estimated,  and  the 
object’s  location  is  approximated  using  the  method  presented  here.  The  estimated  zc  location  of  the 
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Figure  2.25.  Estimated  zc  for  the  8 1-mm  UXO. 
The  true  zr  =  60cm. 
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object’s  center  is  depicted  on  Figure  2.25,  and  the  good  agreement  indicates  that  the  proposed 
approach  is  applicable  for  TD  EMI  data  as  well.  These  results  were  obtained  using  only  the  first 
time  channel,  but  our  studies  indicate  that  other  time  channels  produce  similar  results.  Multiple  data 
channels  can  simultaneously  be  used  to  obtain  more  reliable  estimates  of  the  object’s  position. 

2.8. 7.c.  MPV-TD  data  blind-test  analysis 

We  also  tested  the  NSMS  technique  and  the  HAP  approach  using  TD  MPV  sensor  data. 
The  MPV  TD  sensor  consists  of  two  transmitter  loops  with  3 7. 5 -cm  radii  and  five  triaxial  receiver 
loops/cubes  (see  Figure  26).  The  receivers  are  located  as  follows:  Cube  #0  upper  center;  Cube  #1 
lower  center;  Cube  #2  left  of  center  (x  =  -39.6  cm);  Cube  #3  forward  of  center  (y  =  39.6  cm);  and 
Cube  #4  right  of  center  (x  =  39.6  cm).  These  receivers  accurately  measure  the  complete  transient 
response  over  a  wide  dynamic  range  of  time  going  from  100  ps  to  25  ms.  The  measurements  were 
conducted  by  G&G  personnel  at  the  UXO  test  stand  of  the  USACE-ERDC  Laboratory  in 
Vicksburg,  MS.  The  measurement  platform,  which  is  made  of  nonmetallic  fiberglass  materials,  has 
a  usable  measurement  area  of  3  x  4  m  in  extent.  The  sensor  is  mounted  on  a  robotic  arm  that  can  be 
moved  and  controlled  around  the  measurement  area  using  software  running  on  a  personal  computer. 
Automated  remote  controls  are  used  to  position  ordnance  items  at  an  accurate  depth  (to  within 
1  cm)  below  the  measurement  area.  The  sensor  can  be  positioned  with  an  accuracy  of 
approximately  1  mm.  All  data  presented  here  were  collected  by  personnel  from  G&G  Sciences  Inc. 
The  data  were  collected  on  an  89  point  grid  for  objects  at  different  orientations  and  depths.  The 
response  of  each  object  was  represented  with  only  five  NSMS  belts. 


Figure  2.26.  Actual  and  schematic  diagram  of  the  MPV-TD  sensor. 


32 


Time  channel  #1 


Figure  2.27.  Predicted  and  actual  TD  MPV  data  for  a  60-mm  UXO. 


We  first  test  the  accuracy  of  the  NSMS  technique  against  well  controlled  test-stand  data  for 
a  60-mm  UXO.  Figure  2.27  compares  predicted  and  actual  data  at  a  fixed  time  gate  for  the  x 
components  of  the  scattered  magnetic  fields  as  measured  by  three  receivers  (Cubes  #0,  1,  and  2;  see 
Figure  2.26).  The  results  clearly  show  that  the  NSMS  can  predict  the  EMI  response  of  the  60-mm 
UXO. 


Once  the  accuracy  of  the  NSMS  has  been  established,  we  proceed  to  apply  the  technique  to 
TD-MPV  blind  test  data.  These  blind  test  data  were  collected  for  five  unidentified  targets  at 
CRREL  facilities  by  CRREL  personnel  using  the  MPV.  The  test  objects  considered  were  five  UXO 
(a  BLU-26  bomblet,  105-mm  and  81 -mm  projectiles,  a  60-mm  mortar,  and  a  57-mm  bullet)  for 
which  very  detailed  data  were  collected  under  controlled  conditions  at  the  USACE-ERDC  test  site 
in  Vicksburg.  Using  those  data  sets  we  determined  the  total  NSMS  Q(t )  values  for  the  known  UXO. 
Target  EMI  responses  were  measured  over  34  grid  points  at  two  elevations,  but  in  this  study  only 
data  at  one  elevation  were  used.  First,  the  locations  and  orientations  of  the  objects  were  determined 
using  the  new  algorithm.  The  true  locations  and  orientations  of  the  five  targets  are  given  in  Table 
2.3;  the  inverted  estimates  are  summarized  in  Table  2.4. 

The  comparison  between  inverted  and  true  values  for  position  and  orientation  shows  that 
the  predicted  values  in  general  are  very  close  to  the  true  values.  Sometimes  the  azimuthal  angle  (j)  is 
seen  to  be  flipped;  this  is  mainly  due  to  BOR  symmetry  and  does  not  constitute  an  error. 

Once  the  targets’  locations  and  orientations  were  determined  we  computed  the  total  NSMS 
for  the  five  objects  and  compared  them  to  the  five  total  NSMS  values  stored  in  the  library  of  UXO 
items.  Comparisons  between  the  total  NSMC  Target  #1  and  those  of  the  library  UXO  are  depicted 
in  Figure  2.28.  Note  that  the  NSMS  library  was  created  using  TD  MPV  data  collected  at  the  UXO 
test-stand  site  when  the  objects  were  oriented  horizontally  (dip  =  90°)  with  respect  to  the 
measurement  grid.  The  comparison  shows  that  the  total  NSMS  (a  scalar  discriminator)  is  directly 
related  to  the  object’s  size  and  material  properties:  as  the  size  of  the  object  increases  so  does  the 
total  NSMS  amplitude.  In  addition,  the  inverted  total  NSMS  for  Target  #1  correctly  identified  the 
actual  target.  This  demonstrates  that  the  total  NSMS  is  independent  of  the  target’s  location  and 
orientation  and  is  characteristic  of  the  object.  In  this  test  we  used  TD-MPV  data  collected  at  17 
points  and  one  elevation. 
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Table  2.4.  Ground  truth  for  the  MPV-TD  blind-test  data  runs. 


Target 

ID 

xO  (cm) 

yO  (cm) 

zO  (cm) 

<K°) 

dip  ( ) 

i 

81-mm 

-23.26 

22.5 

56.16 

0 

18.3 

2 

105-mm 

-20.26 

22.5 

69.14 

180 

18.6 

3 

BLU-26 

0.00 

22.5 

43.21 

0 

0.0 

4 

57-mm 

5.22 

22.5 

51.45 

180 

306.6 

5 

60  -mm 

0.00 

22.5 

54.50 

0 

0.0 

Table  2.5.  Inverted  results  for  position  and  orientation  for  MPV-TD  blind-test  data. 
Numbers  with  an  asterisk  are  arbitrary  due  to  BOR  considerations. 


Target 

ID 

xO  (cm) 

y0  (cm) 

zO  (cm) 

d>(°) 

dip  (°) 

1 

81-mm 

-21.16 

23.6 

56.56 

0 

-2 

2 

105  mm 

-5.98 

23.7 

67.57 

170 

-20.8 

3 

BLU-26 

0.25 

22.8 

47.55 

*0 

*180.0 

4 

57-mm 

0.75 

22.6 

58.19 

0 

320.4 

5 

60-  mm 

0.67 

19.3 

54.02 

*3 

0.0 

Figure  2.28.  NSMS  model  comparison  to  library  for  Target  #1. 
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Figure  2.29.  NSMS  model  comparison  to  library  for  Target  #4  (left)  and  Target  #5  (right). 


The  algorithm  identified  correctly  all  remaining  targets.  It  is  worth  noticing,  in  the  case  of 
Targets  #4  and  #5,  that  the  total  NSMS  was  able  to  differentiate  between  the  57-mm  and  60-mm 
UXO  even  though  they  have  comparable  geometrical  sizes  (Figure  2.29).  The  results  show  that  at 
early  times/high  frequencies  the  total  NSMS  from  both  objects  are  similar,  due  to  the  small  skin 
depth;  however,  at  late  times/low  frequencies  the  skin  depth  becomes  significant  and  the  total 
NSMS  of  sizable  objects  is  separable  and  depends  on  the  object’s  material  properties  and  metal 
content,  demonstrating  a  good  discrimination  capability.  To  further  illustrate  how  the  object’s  EMI 
response  affects  the  total  NSMS,  on  Figure  2.30  we  added  the  total  NSMS  of  a  non-permeable 
sphere.  The  non-permeable  sphere  has  the  same  size  as  the  BLU-26,  which  is  made  of  permeable 
steel.  The  results  show  that  at  early  times  the  NSMS  amplitudes  of  the  two  spherical  targets  are 
similar;  at  late  times,  two  objects  of  similar  geometric  size  but  different  material  properties  have 
different  total  NSMS. 
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Figure  2.30.  NSMS  model  comparison  to  library  for  Target  #2  (left)  and  Target  #3  (right). 


Finally  we  investigate  the  data  density  requirement  issue  for  UXO  discrimination.  UXO 
cleanup  demands  real/near-real  time  detection  and  discrimination  systems  to  keep  cleanup  costs 
reasonably  low  by  minimizing  the  time  spent  on  each  anomaly;  moreover,  such  systems  should 
discriminate  between  UXO  and  non-UXO  items  accurately  and  reliably.  To  address  this  problem 
here  we  study  the  TD-MPV  data-density  requirements  for  the  NSMS  model.  Studies  were  done  for 
Targets  #2  and  #3.  The  targets’  total  NSMS,  as  well  as  their  positions  and  orientations,  were 
inverted  using  either  17  or  5  measurement  points.  The  inverted  total  NSMS  values  as  functions  of 
time  for  both  cases  are  similar  and  they  correctly  coincide  which  the  true  object’s  total  NSMS.  Thus 
the  results  suggest  that  the  amplitude  of  the  total  NSMC  remains  stable  as  the  number  and  density 
of  data  points  are  reduced  substantially  for  single-object  discrimination. 
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2.8.7.d.  EMI  Array 


Finally  the  HAP  method  was  adapted  to  the  recently  developed  time  domain  EMI  sensor 
array.  The  sensor  consists  of  25  transmit/receive  pairs  in  a  5  x  5  array  of  time-domain  EMI  sensors 
with  a  35-cm  diameter  transmitter  loop  and  a  25-cm  receiver  loop.  The  sensor  transmits  signal  at  a 
time  in  sequential  order  for  1  to  25,  and  for  each  transmitter  all  receivers  receive.  The  receivers 
measure  the  complete  transient  response  over  a  wide  dynamic  range  of  time  going  approximately 
from  100  jlis  to  25  ms.  The  sensor  provides  625  spatial  data  with  unprecedented  position  accuracy. 
The  data  were  collected  by  NRL  personnel  at  Blossom  Point.  The  actual  and  schematic  diagram  of 
the  NRL  EMI  sensor  are  shown  on  Figures  2.31  and  2.32  respectively.  The  measured  data 
distribution  for  the  Camp  Sibert  4.2-inch  mortar  are  illustrated  on  Figures  2.33  and  2.34.  On  these 
figures  each  subplot  corresponds  to  a  fixed  transmitter  in  subsequent  order  from  1  to  25.  The  HAP 
method  was  applied  to  the  Camp  Sibert  4.2-inch  UXO  TD  EMI  array  data.  The  target  was  placed 
under  the  sensor  at  x0=  y0  =0,  z0  =  -60  cm.  The  TD  data  were  for  each  transmitter  and  123  time 
gates.  The  object’s  position  was  inverted  for  each  time  channel  and  transmitter  using  the  HAP 
technique.  Figure  2.35  shows  inverted  location  versus  transmitter  number  for  each  time  channel. 
The  results  demonstrate  the  superior  performance  of  the  HAP  method  for  predicting  an  object’s 
position. 


Figure  2.31.  EMI  TD  array  sensor  in  operation. 
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Figure  2.32.  Schematic  diagram  of  the  EMI  TD  array  sensor. 
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Figure  2.33.  Measured  EM  field  distribution  on  5  x  5  array  receivers  for  time  channel  #1  for 
different  (from  0  to  24)  transmission;  each  subplot  corresponds  to  a  fixed  trasmitter. 
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Figure  2.34.  Measured  EM  field  distribution  on  5x5  array  receivers  for  time  channel  #40  for 
different  (from  0  to  24  transmission);  each  subplot  corresponds  to  a  fixed  transmitter. 
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Figure  2.35.  Inverted  position  for  Camp  Sibert  4.2-inch  UXO  vs.  transmission  number  at  different 
time  channels. 
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2.9.  Conclusion 


A  new  physics-based  analytic  expression  was  developed  for  estimating  a  buried  object’s 
location,  orientation,  and  magnetic  polarizability  without  solving  the  traditional  ill-posed  inverse 
problem.  The  approach  uses  only  three  global  quantities:  the  magnetic  field  H,  the  vector  potential 
A,  and  the  scalar  potential  y/.  The  algorithm  was  tested  using  both  time-  and  frequency-domain 
synthetic  data  as  well  as  actual  EMI  data  for  dipoles,  spheres  and  UXO,  with  and  without  added 
noise.  The  approach  is  valid  for  both  land-based  and  underwater  UXO  scenarios.  The  2D 
normalized  surface  magnetic  source  (NSMS)  model  has  been  developed  to  find  an  estimate  of  the 
scattered  magnetic  field  and  vector  and  scalar  potentials  from  magnetic  field  data  in  the  case  of 
monostatic  data.  The  method  is  a  generalized  surface  dipole  model  approach  for  approximating 
EMI  responses  from  highly  conducting  and  permeable  metallic  objects.  The  procedure  of 
determining  the  NSMS  amplitudes  distributed  on  a  fictitious  surface  just  beneath  the  measurement 
grid  points  is  realized  in  a  collocation  based  algorithm  from  monostatic  as  well  as  multistatic  data. 

Blind-test  analyses  were  carried  out  for  the  MPV-TD  and  GEM-3D+  sensors  using  the 
normalized  total  surface  magnetic  source  model  in  order  to  classify  between  UXO  and  non-UXO 
items.  In  these  tests  the  NSMC  model  was  combined  with  a  physics-based  approach  for  estimating 
buried  objects’  positions  and  orientations.  During  the  blind-test  UXO  discrimination  study  we  first 
estimated  object  positions  and  orientations  via  the  HAP  technique.  Then  we  determined  the 
amplitudes  of  the  total  NSMS  from  data  for  each  unknown  target,  and  finally  the  inverted  total 
NSMC  time/frequency  dependent  curves  were  used  as  discriminators  by  comparing  them  with  the 
total  NSMS  of  known  UXO  types  from  a  pre-computed  library.  In  all  cases  reported  here  the 
method  was  able  to  identify  with  reasonable  accuracy  the  targets’  positions  and  orientations,  and  at 
the  end  the  NSMS  was  able  to  find  the  correct  UXO  among  the  possible  ordnance  and  clutter  items. 

The  HAP  methods  were  extended  for  the  NRL  EMI  array  TD  sensor  as  well.  The  results 
reported  here  clearly  demonstrate  the  superior  performance  of  the  HAP  technique  for  estimating 
location  and  orientation  of  single  or  aggregate  buried  targets  when  using  actual  EMI  data  from  the 
NRL  TD  EMI  array,  MPV-TD,  GEM-3 D,  and  EM63  instruments.  Moreover,  the  method  is  robust 
with  respect  to  noise. 
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Chapter  3 


Pole-series  and  pseudospectral  approaches  for 
estimating  a  buried  object’s  depth 


3.1.  Introduction 

In  this  chapter  we  consider  the  pole-series-expansion  (PSE)  and  pseudospectral 
approaches  to  estimate  the  locations  and  orientations  of  buried  objects.  The  PSE  uses  measured 
scattered  electromagnetic  (EM)  field  information  at  a  given  set  of  points  to  provide  a  very  fast 
approximation  for  estimating  the  scattered  EM  field  back  to  its  origin  at  any  spatial  resolution. 
The  method  requires  a  bi-static  scattered-field  distribution  on  the  spatial  grid.  Since  most  if  not  all 
sensors  output  mono-static  scattered  fields,  here  we  combine  the  pole-series  approach  with  the 
normalized  surface  magnetic  source  (NSMS)  method  to  predict  buried  objects’  locations  and 
orientations:  initially  the  NSMS  is  used  to  generate  bi-static  EMI  data  from  measured  mono-static 
data,  and  then  the  pole-series  expansion  is  employed  to  localize  scattered  field  singularities 
(SFS),  i.e.,  to  find  the  objects. 

Another  technique  considered  here  for  locating  buried  objects  is  the  pseudospectral  finite- 
difference  (PSFD)  method,  which  is  derived  from  spatial  finite-difference  (FD)  approximations. 
In  the  magnetoquasistatic  regime,  this  technique  is  based  on  a  volumetric  sampling  of  the 
magnetic  field  within  and  around  the  region  of  interest.  In  the  standard  FD  technique,  the  spatial 
sampling  is  set  by  the  user  at  a  sub-spatial-wavelength  resolution  in  order  to  sample  properly  the 
highest  spatial  frequencies  thought  to  be  important  in  the  physics  of  the  problem.  This 
requirement  on  the  space  discretization  has  allowed  the  successful  application  of  FD  methods  to 
solve  Laplace’s  equation  for  a  wide  range  of  electromagnetic  problems  using  second-order- 
accurate  central  differences  on  a  staggered,  uncollocated  Cartesian  lattice.  To  determine  the 
location  of  the  SFS  in  the  EMI  frequency  regime  it  is  necessary  to  solve  simultaneously  the  curl- 
free  and  divergence-free  Maxwell  equations  for  the  scattered  magnetic  field,  which  requires  the 
use  of  different  staggered,  uncollocated  Cartesian  space  lattices.  This  Yee-type  space  stepping  is 
the  source  of  numerical  dispersion  errors.  To  minimize  the  dispersion  error  in  the  computation  of 
spatial  derivatives  we  employ  the  pseudospectral  FD  method  along  the  orthogonal  x-  and 
y-directions.  The  pseudospectral  time-domain  (PSTD)  method  calculates  spatial  derivatives  using 
the  differentiation  theorem  for  Fourier  transforms;  this  process  converges  with  infinite-order 
accuracy  for  grid-sampling  densities  of  two  or  more  points  per  spatial  wavelength.  Our  procedure 
models  magnetic-field  derivatives  along  the  x-  and  y-directions  with  the  PSTD  method  while 
keeping  second-order  central  differences  on  a  staggered,  uncollocated  Cartesian  lattice  along  the 
z-axis;  this  allows  the  refocusing  of  the  scattered  magnetic  field  backward  in  space. 


3.2.  The  pole-series-expansion  approach  in  EMI 


This  method  for  estimating  the  depth  of  a  buried  object  requires  an  essentially  continuous 
representation  of  the  object’s  response  over  a  specific  spatial  range  above  the  measurement 
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surface.  Measuring  the  scattered  EM  field  with  sufficient  resolution  can  be  very  expensive  and 
time-consuming,  especially  when  performing  real-field  UXO  detection  and  discrimination. 
However,  this  problem  can  be  overcome  by  employing  the  NSMS  technique. 

The  spatial  distribution  of  a  linear  system’s  EMI  response  may  be  optimally  represented 
by  Cauchy’s  method.  Consider  the  rational  expansion 


F(s)  =  -^-1  +  Error(s)  = 
D(s) 


+  Error(s) , 


(3.1) 


where  F  is  the  response  of  the  system  and  5  the  coordinate.  The  error  in  this  approximation 
depends  on  the  maximum  orders  n  and  d  of  the  power  series  expansions  in  the  numerator  N  and 
denominator  D  as  well  as  on  the  method  used  to  determine  the  coefficients  Nt  and  Dt.  A  simple 
technique  to  compute  these  parameters  is  to  multiply  (3.1)  by  D  and  evaluate  it  at  a  set  of  m  depth 
points  sk : 


F(sk)EDisL-ENiS'k  =  Error(sO£DisL  =  Ek,  k  =  l,...,m,  (3.2) 

i=0  i=0  i=0 

where  E  is  a  vector  of  unknown  errors.  When  F  is  known  atm>?2  +  <i+l  points  (3.2)  is  a 
linear  system  of  m  equations.  One  then  can  evaluate  the  parameters  D{  and  Nt  in  such  a  way  that 
the  square  norm  of  the  error  vector  is  minimized.  Before  this  is  done,  one  should  note  that  not  all 
parameters  are  independent  because  the  numerator  and  denominator  in  (3.1)  may  be  scaled  by  an 
arbitrary  factor.  For  this  reason,  one  of  the  parameters  may  be  set  equal  to  1.  It  is  reasonable  to  set 
Dd=  1 .  One  then  obtains 


F(sk)XD,sL  "SNiSi  =  -F(sk)sJ+E„=Rk+Ek>  k  =  l,...,m,  (3.3) 

i=0  i=0 

where  R  is  a  known  right-hand-side  vector.  Note  that  (3.2)  can  be  solved  in  such  a  way  that  the 
error  vector  E  is  zero  when  m=Mo  =  n  +  d+  1,  because  one  then  obtains  a  square  matrix  system. 
This  does  not  imply  that  Error(s)  becomes  zero  as  well.  Especially  when  the  sample  values  F(sk ) 
are  only  approximately  known — which  is  always  the  case  in  practice — it  is  more  reasonable  to 
work  with  an  overdetermined  system  of  equations  in  which  m  >  m0.  Reasonable 
overdetermination  should  implicitly  provide  a  “noise”  minimization.  The  most  difficult  problem 
is  to  determine  the  required  maximum  orders  n  and  d  of  the  power  series  of  the  numerator  and 
denominator.  Both  depend  very  much  on  the  size  of  the  space  of  interest,  the  desired  accuracy, 
and  the  complexity  of  the  system.  Since  the  scattered  field  has  l/R3  singularities,  it  is  reasonable 
to  limit  the  maximum  orders  n  and  d  by  a  value  typically  not  higher  than  6  and  to  subdivide  the 
spatial  interval  into  two  or  more  additional  parts  when  the  PSE  approximation  is  not  accurate 
enough.  The  algorithmic  block  scheme  of  the  PSE  procedure  used  for  predicting  object  depths  is 
presented  on  0.  The  method  is  adaptive  and  starts  with  small  orders — -i.e.,  n  and  d  values — and 
with  a  small  number  of  test  points  according  to  an  overdetermination  factor  specified  by  the  user. 
It  then  increases  the  order  by  1  and  compares  the  resulting  PSE  approximations.  When  the 
differences  between  the  two  approximations  are  below  a  user-defined  error  bound  over  the  entire 
burial  depth  range  and  when  all  s  parameters  are  within  the  range  [0,1],  the  PSE  approximation  is 
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good  enough  and  the  procedure  stops.  Otherwise,  the  algorithm  inserts  new  test  points  within  the 
depth  range  as  follows:  if  an  s  parameter  is  out  of  the  range,  it  inserts  a  new  test  point  at  the  depth 
where  the  biggest  distance  from  the  range  [0,1]  is  encountered.  Otherwise,  it  searches  for  the 
maximum  difference  between  the  current  and  previous  PSE  approximation  and  sets  the  new  test 
point  at  the  corresponding  depth. 


Figure  3.1.  Schematic  diagram  of  the  PSE  algorithm. 
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3.3.  Pseudospectral  FD  scheme  in  EMI 


The  objective  of  this  section  is  to  outline  the  basic  theory  of  the  pseudospectral  finite- 
difference  (PSFD)  method  to  invert  for  the  location  of  a  buried  object  via  a  finite-difference  (FD) 
backward  scheme.  The  first  step  of  the  proposed  technique  is  to  obtain  measured  EMI  data  at  a 
plane  with  z  =  h,  x  e  [- a ,  a\,  y  e  [- b ,  b\.  Then  the  entire  volume  is  divided  into  subvolumes/cells 
characterized  by  the  triple  index  (7,  j,  k)  with  I  =  1,  . . .,  Nx,j  =  1,  . . Ny,  k  =  1  ,...,NZ.  The  PSFD 
technique  is  combined  with  the  NSMS.  First  the  NSMS  amplitudes  are  determined  from  the 
measured  data  at  the  plane  zNz  =  h ,  xz  =  -a  +  dx  (I-  1),  yj  =  -b  +  dy(j  -  1)  (where  dx,  dy ,  and  dz 
give  the  size  of  the  (/,  j ,  A:)  cell)  by  solving  a  linear  system  of  equations.  At  the  next  step  a  scalar 
magnetic  field  potential  \| /  is  generated  via  SMS  on  the  plane  zNz  =  h  +  dzl 2,  xz  +  dxl 2,  yy  +  dyl 2 
(/=  1  j  =  1,...,AV.).  We  can  then  use  the  scalar  magnetic  potential  \|/  at  zNz  =  h  +  dzl 2, 

X/  +  dxl 2,  +  dyl 2  and  the  magnetic  field  //z  at  z/Vz  =  /z,  xz  +  dxl 2,  +  dyl 2  to  recover  both  the 
magnetic  field  and  the  scalar  potential  by  solving  H  =  -V\|/  and  V-H  =  0  at  each  (/,  j,  k)  cell. 

To  illustrate  the  technique,  let  us  consider  an  object  placed  in  free  space  and  illuminated 
with  a  primary  magnetic  field.  The  primary  magnetic  field  penetrates  the  object  and  induces  eddy 
currents  within  it;  these  eddy  currents  in  turn  produce  a  measurable  secondary  magnetic  field 
outside  the  object.  Our  goal  is  to  use  the  magnetic  fields  at  two  surfaces  (x-y,  z  =  h  and  x-y, 
z  =  h  +  dzl 2)  to  predict  the  scattered  magnetic  field  below  the  surface  x-y,  z  =  h.  To  do  that  let  us 
combine  the  x  and  y  components  of  the  curl-free  Maxwell  equation  for  the  scattered  magnetic 
field  Hsc  with  the  divergence-free  equation: 


Figure  3.2.  Yee  FD  scheme  for  the  PSFD  approach. 
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dHsxc  _  dHszc 
dz  dx 
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m; 

dz 


dHs; 

dz 


(dHsc  dHs;) 

— —  + — — 

8x  dy  J 


(3.4) 


By  using  the  standard  FD  scheme  (3.4),  the  curl-  and  divergence-free  Maxwell  equations  for  Hsc 
in  FD  form  can  be  written  as 


H?  (i  +  0.5,  j,k)~  H*  U  +  0.5,  j,  k  - 1) 
Az 

HI  (i,  j  +  0.5,  k)  -  HI (i,  j  +  0.5,  k  - 1) 
Az 


0-5)  -  H^_  (/  - 1,  j,  k  +  0.5) 
Ax 

HI  (i,  j,  k  +  0,5)  -  H*  (j,  j  -l,  k  +  0.5) 
Ay 


(a) 

(b) 


H^  (J  +  0-5,  j  +  0.5,  k)  -H“  (i  +  0.5,y  +  0-5,  k  -1)_ 

Az 

f  Hsxc(i,j  +  0.5, k  +  0.5) - Hsxc(i -1J  +  0,5, k  +  0.5)  |  H*(i +  0.5,  j,k +  0.5)  - Hsyc(i  +  0.5,  j  -l,k  +  0.5) A 


V 


Ax 


Ay 


(C) 


(3.5) 


where  the  subscripts  /,  /,  and  k  are  the  Cartesian  coordinate  indices  of  the  spatial  sampling  point 
(I  Ax,  j  Ay,  k  Az)  in  the  three-dimensional  (3-D)  lattice.  According  to  the  standard  FD  scheme  the 
curl-free  equations  (3.5)(a)  and  (3.5)(b)  and  the  divergence-free  equation  (3.5)(c)  require 
difference  lattice  points  (see  Figure  3.2:  black  and  red  arrows  apply  to  the  curl  and  divergence 
equations  respectively).  This  complicates  the  applicability  of  the  standard  FD  scheme  to  refocus 
the  scattered  magnetic  field  backwards  to  its  origin.  To  overcome  this  difficulty  and  avoid  the 
spatial  derivatives  on  a  staggered,  un-collocated  grid,  we  employ  a  pseudospectral  FD  (PSTD) 
technique  for  calculating  the  derivatives  with  respect  to  the  x-  and  y-coordinates.  The  PSTD 
method  utilizes  discrete  Fourier  transforms  to  evaluate  spatial  derivatives  on  unstaggered, 
collocated  grids.  For  example,  the  x-derivative  of  a  general  field  component  that  is  known  at 
all  grid  points  /  along  the  x-direction  can  be  computed  as  follows: 


d' P 


dx 


(3.6) 


where  kx  is  the  Fourier-transform  variable  representing  the  x-component  of  the  numerical  wave 

vector  and  F  and  F-1  denote,  respectively,  the  forward  and  the  inverse  discrete  Fourier 
transforms  along  the  x-direction.  The  Fast  Fourier  Transform  (FFT)  algorithm  is  used  to 
implement  the  approximation.  The  potential  limitation  caused  by  the  periodic  boundary 
conditions  inherent  in  the  FFT  can  be  eliminated  by  use  of  absorbing  boundary  conditions.  Note 
that,  according  to  the  Nyquist  sampling  theorem,  the  representation  of  approximation  (3.6)  is 
exact  for  grid-sampling  densities  of  two  or  more  cells  per  wavelength.  Hence  any  dispersion 
errors  in  the  PSFD  method  introduced  here  arise  only  from  the  second-order-accurate  central 
differences  on  the  staggered,  uncollocated  Cartesian  space  lattice  along  the  z-direction.  We  note 
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that  there  is  a  potential  problem  with  the  PSFD  for  focused,  compact,  delta-function  type  sources, 
since  this  type  of  source,  which  represents  a  spatially  discretized  delta  function,  may  cause 
difficulties  when  using  the  FFT.  However,  we  can  avoid  this  difficulty  by  employing  the  pole- 
series  expansion  technique  described  above.  After  combining  (3.5)  and  (3.6),  the  magnetic  field 
at  height  z^_i  can  be  expressed  with  the  magnetic  field  and  its  derivative  on  the  unstaggered  grid 
as 


Hx  (i,  j,  k  -  0.5)  =  Hx  (i,  j,  k  +  0.5)  -  Az  •  F-1  [- jkxF(Hz )] 

Hy  (i,  j,  k  -  0.5)  =  Hy  (i,  j,  k  +  0.5)  -  Az  •  F1  [- jkyF(Hz )]  .(3.7) 

Hz  (i,  j,  k  -  0.5)  =  Hz  (i,  j,  k  +  0.5)  +  Az  •  [F-1  [-jkxF(Hx )]  +  F"1  [- jkyF(Hy  )]] 


3.4.  Noise  analysis  in  the  PSE  and  PSFD  approaches 

The  relation  between  a  measured  data  point  dj  and  the  prediction  of  a  forward  model  at 
the  same  point  (which  we  denote  Fj)  can  be  expressed  as 

dj  =Fj(m)  +  8j,  (3.8) 

where  j  is  one  of  N  data  points,  m  is  a  vector  of  model  parameters,  and  Sj  is  the  error  of  the  j- th 
datum.  In  general  the  error  can  fall  into  one  of  two  categories:  a)  those  related  to  the  uncertainties 
in  modeling  parameters  (like  sensor  positioning  and  orientation)  not  included  in  the  model  vector 
m  which  produce  errors  in  the  predicted  data,  and  b)  those  that  are  strictly  related  to  the  sensor, 
even  in  the  absence  of  a  target:  the  effects  of  instrument  drift  and  electronics  noise,  among  others. 
All  of  these  errors  in  the  data  enter  (3.8)  as  an  additive  term,  and  usually  appear  in  the  inverse 
problem  through  the  data  covariance  matrix.  The  data  covariance  matrix  adjusts  dj ,  the  relative 
contribution  of  each  data  point  to  the  objective  function,  and  therefore  controls  how  closely  each 
datum  is  fit  by  the  predicted  data.  In  our  subsequent  analysis  we  assume  independently 
distributed  Gaussian  random  errors. 


3.5.  Results 


3.5.1  PSE  applied  to  a  dipole 

We  designed  various  studies  to  demonstrate  the  applicability  of  the  PSE  approach  for 
determining  SFS  and  hence  for  estimating  the  location  of  a  buried  object.  We  initially  placed  a 
magnetic  dipole  at  the  origin  of  a  Cartesian  coordinate  system  and  calculated  its  magnetic  field 
from  30  cm  to  40  cm.  We  used  the  field  values  in  this  segment  to  determine  the  PSE  expansion 
coefficients  and  then  predict  the  magnetic  field  at  a  segment  between  -60  cm  and  60  cm.  A 
comparison  between  the  PSE  predictions  and  the  actual  fields,  depicted  on  Figure  3.3,  shows  that 
the  PSE  has  a  reasonable  predicting  power. 
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Magnetic  dipole  radiation 


Magnetic  dipole  radiation 


Figure  3.3.  Comparison  between  predicted  and  actual  fields  for  a  magnetic  dipole. 
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Distance  [cm] 


Distance  [cm]  Distance  [cm] 


Figure  3.4.  (a)  Experimental  setup  for  Hz  field  upper  continuation  and  (b)  clutter  items  used  in 
the  experiment.  Plots:  Inphase  (top)  and  quadrature  (bottom)  Hz  field  distribution  on 
the  measurement  plane  as  predicted  by  NSMS  (left)  and  measured  by  the  GEM-3 
(right). 
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3.5.2.  Extending  the  EMI  field  upwards 


In  this  section  we  present  results  of  further  numerical  experiments  carried  out  to  assess 
the  applicability  of  the  combined  NSMS/PSFD  method  for  determining  the  location  of  a  buried 
object  without  having  to  solve  a  time-consuming  minimization  problem.  First  we  compare  an 
extended  EMI  magnetic  field  computed  via  NSMS  to  measured  frequency-domain  data.  In  these 
tests,  the  NSMS  model  was  used  to  extend  monostatic  EMI  data  above  the  measurement  plane 
(which  corresponds  to  the  bottom  plane  of  the  “computational  space”  shown  on  Figure  3.4).  The 
scattered  magnetic  field  is  produced  by  a  set  of  magnetic  charges  placed  on  a  fictitious  surface. 
This  fictitious  surface  is  divided  into  small  patches  (Figure  3.4),  and  at  each  subsurface  the 
normal  component  of  the  primary  magnetic  field  is  determined.  The  amplitudes  of  the  magnetic 
charges  at  each  subsurface  are  scaled  by  the  normal  component  of  the  actual  primary  magnetic 
field.  The  NSMS  amplitudes  are  determined  by  matching  the  modeled  magnetic  field  to  the 
measured  field  at  a  selected  set  of  grid  points  on  the  measurement  plane.  Once  these  amplitudes 
are  found,  the  scattered  EMI  field  is  extended  into  the  computational  space.  Note  that  the  NSMS 
sources  are  distributed  on  a  fictitious  plane  surface  placed  between  the  ground  and  the 
computational  space,  corresponding  to  the  red  grid  on  0.  In  this  arrangement  the  NSMS  reproduce 
the  magnetic  fields  radiated  by  all  the  metallic  objects  distributed  beneath  the  fictitious  surface. 
For  NSMS  validation,  we  distributed  various  clutter  items  (shown  in  0)  on  the  plane  surface  and 
collected  data  at  7  x  7  grid  points  on  two  surfaces  using  a  GEM-3  sensor.  The  GEM-3  consists  of 
two  parallel  concentric  circular  transmitter  loops  with  currents  flowing  in  opposite  directions  and 
scaled  so  that  the  total  primary  field  vanishes  at  the  common  center,  where  a  receiver  coil  is 
located.  The  transmitter  loops  are  idealized  as  infinitely  thin  line  sources  with  radii  ax  and  a2  and 
currents  7j  and  I2  =  -I\  a2lax.  The  NSMS  were  distributed  on  a  plane  surface  5  cm  above  the 
clutter  surface.  Matching  the  NSMS-modeled  magnetic  field  to  the  actual  data  at  h  =  10  cm  above 
the  clutter  surface  enabled  us  to  determine  the  amplitudes  of  the  NSMS,  which  we  then  used  to 
compute  the  scattered  mono-static  field  at  7  x  7  points  on  a  plane  15  cm  above  the  clutter  surface. 
The  GEM-3  data  were  collected  on  the  same  surface  and  using  the  same  grid  points.  The  actual 
GEM-3  data,  at  210  Hz,  are  depicted  on  the  right  column  of  Figure  3.4;  the  left  column  shows  the 
z-component  of  the  magnetic  field  as  predicted  by  the  NSMS  model.  We  find  very  good 
agreement  between  the  predicted  field  and  the  measured  data,  which  demonstrates  that  the  NSMS 
model  can  be  used  successfully  to  extend  EMI  data  in  a  given  computational  space. 


3.5.3.  Combining  NSMS  and  PSE 

A  nice  feature  of  the  NSMS  technique  is  that  it  can  be  used  to  generate  both  monostatic 
and  bistatic  data  for  any  given  sensor.  Bistatic  data  is  important  when  determining  locations  and 
orientations  of  buried  objects  using  the  pole-series  expansion  technique  described  here.  However, 
most  if  not  all  current  state-of-the-art  EMI  sensors  are  monostatic.  The  NSMS  method  readily 
clears  this  hurdle  by  performing  an  accurate  data  conversion  from  monostatic  to  bistatic.  To 
illustrate  this  we  solve  the  complete  EMI  scattering  problem  for  a  highly  permeable  and 
conductive  cylinder  using  a  method  of  auxiliary  sources/thin  skin  approximation  (MAS-TSA) 
code  and  use  the  results  to  determine  the  MAS-based  standardized  excitation  approach  (MAS- 
SEA)  coefficients.  We  then  generate  MAS-SEA  monostatic  data  on  the  gridx  e  [-60  cm,  60  cm], 
ye  [-60m,  60cm],  z  =  45cm  at  10-cm  increments  for  the  cylinder  subject  to  the  primary 
magnetic  field  of  the  GEM-3  sensor.  In  the  next  step  we  distribute  NSMS  sources  on  a  flat 
surface  with  x  e  [-80  cm,  80  cm],  y  e  [-80  cm,  80  cm],  z  =  20  cm,  that  is  conformal  to  but  does 
not  coincide  with  the  measurement  surface.  In  this  arrangement,  the  NSMS  reconstructs  the 
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Distance  [cm] 


magnetic  fields  radiated  by  all  metallic  objects  distributed  in  the  space  z  <  0.  The  amplitudes  of 
the  NSMS  are  determined  by  solving  a  linear  system  of  equations  and  enable  us  to  simulate  the 
bistatic  scattered  EMI  field  on  the  grid,  which  we  can  then  compared  with  the  data  generated  by 
the  MAS-SEA  code.  Figure  3.5  shows  one  such  comparison,  both  on  a  grid  (panels  (a)  and  (b)) 
and  along  a  transect  (panel  (d)).  Panel  (c)  shows  the  difference  between  simulated  and  predicted 
data.  We  see  that  the  NSMS  method  can  reconstruct  bistatic  data  from  monostatic  data  very 
accurately. 
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Figure  3.5  Hz  bistatic  magnetic  field  on  a  grid. 
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3.5.4.  PSE-NSMS  applied  to  measured  EM-63  data 


To  demonstrate  the  applicability  of  the  combined  NSMS-pole  series  expansion  method  we  then 
placed  an  actual  UXO  (a  2. 75 -inch  projectile)  in  the  magnetic  field  of  the  EM-63  sensor.  Data 
were  collected  at  the  UXO  test-stand  site  over  a  set  of  grid  points  on  a  surface  with  x  e  [-1.5  m, 
1.5  m],  y  e  [-1.5  m,  1.5  m],  z  =  60  cm.  NSMS  were  distributed  on  a  flat  surface  with  the  same  x- 
y  range  as  the  measurement  surface  and  z  =  0.  In  this  arrangement,  the  NSMS  reconstruct  the 
magnetic  fields  radiated  by  all  metallic  objects  distributed  in  the  space  z  <  0.  The  NSMS 
amplitudes  were  determined  by  solving  a  linear  system  of  equations  and  used  to  extend  the 
scattered  EMI  field  above  the  measurement  surface  at  z  =  60  cm.  In  order  to  illustrate  the  full 
potential  of  the  PSE-NSMS  approach  for  position  determination  we  also  conducted  a  noise 
analysis:  to  the  NSMS-predicted  Hz(x,y,z )  magnetic  field  signals  we  added  random  noise  as 

Hz  (x,  y,  z)  =  H”sms  (x,  y,  z)  +  e(z) ,  (3.9) 

where  H"sms(x,y,z)  is  the  scattered  magnetic  field  above  the  measurement  surface  and  e(z)  is  a 

Gaussian  random  variable.  The  magnetic  field  predicted  by  the  NSMS  model  at  the  first  time 
channel  is  depicted  in  Figure  3.6  for  different  noise  levels  over  the  x-z  plane  with  y  =  0.  The 
scattered  magnetic  field  is  computed  on  a  plane  with  x  e  [-50  cm,  50  cm]  and  z  e  [-50  cm, 
110  cm];  the  target  is  oriented  horizontally.  For  each  fixed  x  we  use  the  field  values  along  the  z- 
axis  to  determine  the  pole  expansion  coefficients  (with  n  =  6,  d  =  6)  and  reconstruct  the  magnetic 
field.  The  computation  is  straightforward  and  very  fast.  The  maximum  of  the  reconstructed  field 
distribution  (i.e.,  the  SFS)  is  found  to  be  near  the  object  for  very  low  noise  levels  (5%  or  less); 
however,  as  the  noise  level  increases  the  algorithm  becomes  unstable  and  the  predicted  field 
distribution  scatters  around.  Similar  tests  were  performed  for  the  same  2. 75 -inch  mortar  inclined 
45°.  The  results,  depicted  on  Figure  3.7,  lead  to  the  same  conclusion  as  before:  the  maximum  of 
the  inverted  field  is  distributed  around  the  object  and  is  sensitive  to  the  noise.  We  can  also  see 
that  the  object’s  orientation  is  very  difficult  to  estimate. 
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Figure  3.6.  Distribution  of  predicted  Hz  magnetic  field  intensity  over  the  x-z  plane  for  a  2.75-inch 
rocket  oriented  horizontally  and  interrogated  by  the  GEM-3  sensor.  The  different 
panels  correspond  to  increasing  noise  levels:  (a)  0,  (b)  0.5%,  (c)  1%,  (d)  2%,  (e)  3%, 
(f)  5%,  (g)  10%,  (h)  15%. 
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Figure  3.7.  Distribution  of  predicted  Hz  magnetic  field  intensity  over  the  x-z  plane  for  a  2.75-inch 
rocket  at  45°  inclination  and  interrogated  by  the  GEM-3  sensor.  The  different  panels 
correspond  to  increasing  noise  levels:  (a)  0,  (b)  0.5%,  (c)  1%,  (d)  2%,  (e)  3%,  (f)  5%, 
(g)  10%,  (h)  15%. 
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Figure  3.8.  Predicted  and  actual  scattered  magnetic  field  of  a  sphere  at  different  heights  above 
the  sphere  center.  The  y-scale  is  logarithmic. 


3.5.5.  PSFD  for  refocusing  the  scattered  magnetic  field 

In  this  final  subsection  we  use  the  PSFD  to  determine  the  location  of  a  buried  object  by 
refocusing  the  scattered  magnetic  field  backwards  to  its  origin.  The  PSFD  technique  requires 
bistatic  data,  which  can  be  readily  computed  from  monostatic  data  by  using  the  NSMS  method  as 
described  above.  We  consider  a  sphere  placed  in  the  primary  magnetic  field  of  the  GEM-3  sensor. 
We  generated  synthetic  data  at  a  set  of  grid  points  on  a  surface  with  xje  [-2  m,  2  m]  at  heights 
h=z  =  30  cm  and  h=z  +  Az  in  the  Cartesian  coordinate  system  of  Figure  3.2.  Figure  3.8 
compares  the  actual  magnetic  fields  of  the  sphere  and  the  predictions  based  on  (3.7)  at  different 
heights  when  Az  =  1  cm.  The  predicted  magnetic  field  matches  the  ground  truth  at  h  =  29,  20,  and 
10  cm.  However,  when  the  observation  line  is  just  above  the  sphere  there  are  significant 
differences  which  arise  because  the  secondary  source  is  very  concentrated.  This  situation  causes 
difficulties  for  the  FFTs  in  the  PSFD  approach  but  can  be  overcome  by  using  the  pole-series 
expansion  method  described  in  this  chapter.  Similar  tests  with  Az  =  2  cm  yielded  identical  results, 
indicating  that  the  PSFD  method  converges. 
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3.6  Conclusion 


In  this  chapter  we  combined  the  pole-series  expansion  and  pseudospectral  finite- 
difference  method  with  the  normalized  surface  magnetic  source  approach  to  invert  for  the 
location  of  a  buried  object  without  solving  a  traditional  ill-posed  optimization  problem.  We  first 
used  NSMS  to  extend  an  EMI  scattered  field.  Numerical  tests  confirmed  that,  as  demonstrated  in 
the  previous  chapter,  the  NSMS  method  accurately  predicts  monostatic  EMI  responses.  We  then 
employed  the  PSE  or  PSTD  to  focus  a  scattered  magnetic  field  back  to  its  origin.  PSE-NSMS 
studies  were  done  for  a  dipole  and  a  2. 75 -inch  mortar  under  different  levels  of  noise.  We  found 
that  the  PSE-NSMS  technique  is  potentially  able  to  estimate  the  depth  of  a  buried  object  to  an 
accuracy  of  plus  or  minus  one  characteristic  length,  which  could  be  used  as  the  initial  guess  in  a 
nonlinear  optimization.  On  the  other  hand,  this  technique  is  very  sensitive  to  the  noise  level, 
which  limits  its  usefulness  for  actual  UXO  detection  and  discrimination  problems.  We  also  tested 
the  PSFD-NSMS  approach  for  a  sphere.  Our  results  demonstrate  that  the  PSFD-NSMS  method 
can  reconstruct  the  scattered  magnetic  field  to  an  accuracy  of  plus  or  minus  one  sphere  radius; 
however,  it  is  difficult  to  make  correct  depth  estimates  with  this  algorithm,  which  moreover  is 
computationally  expensive  and  sensitive  to  noise. 
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Chapter  4 

Scatterer  localization  using  a  left-handed-medium  slab 


4.1.  Introduction 

Back  in  1968,  V.  G.  Veselago  theoretically  investigated  materials  with  simultaneously 
negative  permittivity  and  permeability,  henceforth  referred  to  as  left-handed  materials  (LHMs), 
and  pointed  out  some  of  their  electrodynamic  properties,  such  as  the  reversal  of  Snell’s  law,  of 
the  Doppler  effect,  and  of  Cerenkov  radiation  (Veselago,  1968).  Veselago’s  ideas  did  not  receive 
much  attention  at  the  time  and  were  eventually  forgotten  because  of  the  unavailability  of  LHMs. 
However,  in  1996  the  first  artificial  LHM  was  realized  using  periodic  structures  (Pendry  et  al , 
1996),  and  since  then  intense  study  has  been  conducted  on  the  theory,  experimental  behavior,  and 
potential  applications  of  LHMs.  As  predicted  by  Veselago,  a  LHM  slab  can  refocus 
electromagnetic  waves  from  a  point  source.  J.  Pendry  and  collaborators  (1996)  extended 
Veselago’s  ideas  and  further  predicted  that  evanescent  waves  could  be  amplified  through  an 
LHM  slab  and  that  the  source  information  could  be  reconstructed  at  the  perfect  imaging  point 
without  loss  of  amplitude  (Pendry  et  al.,  1996;  1999).  Similar  tests  have  been  done  for  spatially 
distributed  sources  (such  as  Gaussian  beams)  propagating  through  a  slab  of  a  left-handed 
medium  (Ziolkowski,  2003),  and  it  has  been  reported  that  Gaussian  beams  at  normal  incidence 
can  be  focused  in  a  planar  LHM  slab.  Since  then  most  of  the  research  activity  on  LHMs  has 
focused  on  the  development  of  actual  LHM  devices,  such  as  perfect  lenses,  resonators,  etc.  The 
aim  of  this  project  was  to  investigate  the  possibility  of  applying  LHMs  and  their  refocusing 
properties  as  a  virtual  mathematical  tool  to  help  determine  the  location  and  geometric  parameters 
of  a  buried  object  from  measured  data  without  having  to  solve  an  ill-posed  inverse-scattering 
problem.  We  would  like  to  emphasize  that  this  approach  does  not  require  any  hardware,  since  it  is 
based  exclusively  on  the  mathematical  properties  of  LHMs. 


4.2.  Governing  equations  in  the  EM  wave  regime 

According  to  Maxwell’s  equations,  in  a  LH  medium  characterized  by  8  <  0  and  \i  <  0 
the  plane-wave  relations  between  electric  ( E(r)  =  E0e-jk  r )  and  magnetic  ( H(r)  =  H0e-jk  r ) 
fields  take  the  form 

k  x  E  =  -co  |  (i  |  H  (4.1) 


and 


k  x  H  =  co  |  8  |  E ,  (4.2) 

where  a  time  dependence  eJC0t  is  assumed.  As  a  result,  the  electric  field,  the  magnetic  field,  and 
the  wave  vector  k  form  a  left-handed  triad.  Thus  the  wave  vector  k,  and  therefore  also  the  phase 
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Figure  4.1.  Configuration  of  Poynting  vectors  (S)  and  wave  vectors  (k)  in  a  LHM  half-space. 


velocity  v^,  exhibit  a  sign  opposite  to  that  in  a  conventional  right-handed  medium  (Figure  4.1).  At 
the  same  time,  Poynting ’s  theorem  is  still  given  by 

S  =  ExH*,  (4.3) 

which  indicates  that  the  electric  and  magnetic  fields  and  the  Poynting  vector  S  are  still  related 
through  the  right-hand  rule  (see  Figure  4.1).  Thus  the  Poynting  vector  S,  and  therefore  also  the 
group  velocity  vgr,  point  in  the  same  direction  as  the  propagation  of  energy,  as  in  right-handed 
(RH)  media,  thereby  satisfying  the  causality  requirement.  Combining  observations  of  the 
directions  of  k  and  S  demonstrates  that  in  a  LHM  the  phase  and  group  velocities  are  antiparallel 
(/.  e. ,  of  opposite  sign),  and  that  the  wave  fronts  travel  toward  the  source. 

To  better  understand  EM  scattering  and  refocusing  phenomena  in  LH  materials,  let  us 
consider  a  TE-polarized  wave  (electric  field  perpendicular  to  the  plane  of  incidence)  incident  on  a 
LHM  half-space,  as  shown  in  Figure  4.2.  The  reflected  and  transmitted  waves  are  known  to 
satisfy  the  law  of  reflection  and  Snell’s  law: 


0;  -  9i, 

(4.4) 

nt  sin0t  =11;  sinBj, 

(4.5) 

where  the  index  of  refraction  of  each  medium  (/?=  /,  t)  is  given  by  the  expression 
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Figure  4.2.  Interface  between  free  space  and  a  LHM  showing  a  negative  angle  of  refraction. 


n ;  = 


(4.6) 


where  £o  and  //o  are  respectively  the  permittivity  and  the  permeability  of  free  space.  The  reflection 
and  transmission  coefficients  are  given  by 


R  =  7j1cos0L_-rj1cos0L 
r/tcosdi  +  r/icosdt 


and 


(4.7) 


T=  2/7,  cos  6>. 
r/tcosOi  +  rji  cos  0t 

where  the  wave  impedance  in  medium  ft  is 


(4.8) 


(4.9) 


In  most  cases  below  we  consider  a  left-handed  metamaterial  matched  to  free  space.  This  means 
that  sr  <  0  and  |Lir  <  0  so  that  n  <  0,  and  that  r|t  =  r|i .  Consequently,  for  a  LHM  with  n  =  — 1 

one  finds  R  =  0  and  T=  1.  Moreover,  Snell’s  Law  indicates  that  the  transmitted  angle  is  negative 
for  oblique  incidence  (at  any  angle)  on  the  interface.  The  negative  angle  of  refraction  has  the 
effect  of  refocusing  the  electromagnetic  field  back  to  the  origin.  In  other  words,  the  planar  LHM 
behaves  as  a  perfect  lens. 
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Figure  4.3.  Schematic  diagram  of  a  2D  EM  problem. 


4.3.  Boundary  conditions  at  a  free-space/LHM  interface 

Let  us  examine  now  the  boundary  conditions  (BCs)  for  electromagnetic  fields  at  the 
interface  between  free  space  and  a  LH  medium.  The  tangential  components  of  E,  H,  and  k  do  not 
“see”  the  discontinuity  and  are  therefore  transmitted  from  one  medium  to  another  unaffected. 
Mathematically,  the  boundary-value  problem  for  the  tangential  components  can  be  summarized 
as  follows: 


nx(Hpr+Hobj)  =  nxHLHM 
nx(Epr+Eobj)  =  nxELHM 


(4.10) 


where  n  is  a  unit  vector  normal  to  the  surface,  Epr  and  Hpr  are  the  primary  electric  and  magnetic 
fields,  Eobj  and  Hobj  are  the  measurable/recordable  scattered  electric  and  magnetic  fields  due  to  the 
object,  and  elhm  and  hlhm  are  the  total  electric  and  magnetic  fields  inside  the  LHM,  which  we 
wish  to  reconstruct.  In  contrast,  the  normal  components  of  the  electric  and  magnetic  vectors 
undergo  a  change  of  sign  at  the  interface  in  addition  to  the  usual  discontinuity  in  magnitude.  The 
sign  changes  for  En  and  Hn  come  directly  from  the  conditions  S\E\n  =  SjEin  with  s2  <  0  and 
jU\H\n  =  jU2H2n  with  ju2  <  0,  while  the  change  of  sign  for  kn  is  a  consequence  of  the  condition  n  <  0. 
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4.4.  Refocusing  of  the  scattered  EM  field  using  LHM 


Let  us  consider  a  permeable  and  conducting  two-dimensional  metallic  object  buried  in 
soil  and  illuminated  with  a  primary  electromagnetic  field  (Figure  4.3).  (A  time  dependence  eJC0t  is 
assumed  and  suppressed  throughout.)  The  primary  electromagnetic  field  induces  electric  currents 
on  the  scatterer’s  surface,  which  in  turn  produce  a  secondary  or  scattered  field  outside.  The 
scattered  electromagnetic  field  is  measured/recorded  at  a  set  of  measurement  points.  Our  goal  is 
to  use  this  measured  secondary  EM  field  to  determine  the  location  of  the  object.  To  do  so,  we  first 
assume  that  the  measurement  surface  is  a  boundary  between  free  space  and  a  LHM.  Then  the 
scattered  EM  field  is  considered  as  a  primary  EM  field  on  the  virtual/fictitious  free-space/LHM 
boundary.  Third,  the  electromagnetic  boundary-value  problem  is  solved,  and  finally  the  EM  field 
is  mapped  inside  the  LHM.  In  the  regime  considered  here  the  electromagnetic  fields  are  governed 
by  vector  wave  equations.  In  the  2D  case  an  electromagnetic  field  can  be  represented  with  a 
scalar  potential  that  satisfies  the  Helmholtz  equation, 

(V2  +  kp)  U  =  -Jz  5(p  -  p') ,  (4.11) 

where  kp  is  the  wave  number  (measured  in  m"1)  in  medium  / 3  and  U  is  the  scalar  potential.  Thus 
the  problem  is  cast  in  terms  of  solutions  of  the  Helmholtz  equation  inside  each  medium.  This 
boundary-value  problem  can  be  solved  using  the  method  of  auxiliary  sources. 


4.5.  The  method  of  auxiliary  sources  for  LHM  in  the  EM-wave  regime 

The  method  of  auxiliary  sources  (MAS)  is  a  numerical  technique  originally  developed  for 
solving  a  wide  class  of  scattering  and  radiation  electromagnetic  problems  that  recently  has  also 
been  applied  to  various  electromagnetic  induction  (EMI)  problems  (Zaridze  et  al .,  2002; 
Shubitidze  et  al ,  2002;  Anastassiu  et  al ,  2002;  Shubitidze  et  al ,  2004).  It  has  been  widely 
demonstrated  to  be  a  general,  robust,  and  accurate  numerical  method  for  the  study  of  EMI 
scattering  by  highly  conducting  and  permeable  targets,  and  has  been  used  to  solve  complex  large- 
scale  electromagnetic  scattering  and  antenna  problems.  Because  of  its  reduced  computational 
complexity,  the  method  shows  great  potential  for  simulating  realistically  elaborate 
electromagnetic  situations,  particularly  those  involving  multiple  objects.  Briefly,  in  the  MAS, 
boundary  value  problems  are  solved  numerically  by  representing  the  electromagnetic  fields  in 
each  domain  of  the  structure  under  investigation  by  a  finite  linear  combination  of  analytic 
solutions  of  the  relevant  field  equations,  corresponding  to  sources  situated  some  distance  away 
from  the  boundaries  of  each  domain.  The  “auxiliary  sources”  producing  these  analytical  solutions 
are  chosen  to  be  elementary  currents  located  on  fictitious  surfaces  that  usually  conform  to  the 
actual  surface(s)  of  the  structure.  In  practice,  only  points  on  the  auxiliary  and  actual  surfaces  are 
required,  so  it  is  not  necessary  to  resort  to  the  detailed  mesh  structures  required  by  other 
techniques  such  as  the  finite  element  method  and  the  boundary  element  method. 

For  a  simple  one-dimensional  uniform  half-space  with  given  permittivity  and 
permeability  we  require  two  auxiliary  surfaces,  one  inside  and  one  outside  the  scattering  object 
(/.£.,  the  half-  space).  The  fields  outside  the  half-space  (region  1)  are  considered  to  originate  from 
a  set  of  auxiliary  sources  {«/}  placed  inside  the  half-space  (surface  S\mx),  while  the  interior  fields 
(region  2)  arise  from  a  set  of  auxiliary  sources  placed  outside  (surface  S2aux).  The  auxiliary 
currents  {at}  and  {Z?/}  radiate  in  an  unbounded  homogeneous  space  filled  respectively  with  the 
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characteristics  of  region  1  and  region  2.  The  electric  and  magnetic  fields  in  regions  1  and  2  are 
linked  by  the  boundary  conditions  (4.10),  which  are  enforced  at  an  array  of  selected  points  on  the 
physical  surface  S  corresponding  to  the  interface. 

The  electric  field  generated  by  this  current  density  is  described  by 

E (x,z)  =  ~^I0H™ (k0  V(x - x0 )2  +  (z - z0 )2  )  ,  (4.12) 

where  H0(2)  is  the  zeroth-order  Hankel  function  of  the  second  kind  and  £0  and  /Jo  are  the  free-space 
wave  number  and  magnetic  permeability,  respectively.  The  intensity  of  the  radiated  magnetic 
field  is  obtained  via  Faraday’s  law.  The  primary  EM  field  impinges  at  the  interface  between  the 
free  half-space  of  region  1  and  the  dissipative  medium  of  region  2;  consequently,  two  secondary 
EM  fields  are  generated.  The  first  one,  denoted  by  Eobj  and  Hobj,  is  the  scattered  field,  which 
propagates  in  the  free  half-space,  and  the  second  one,  denoted  by  elhm  and  Hlhm,  is  the  refracted 
field  that  propagates  in  the  LH  medium.  The  incident  field  (and  thus  the  entire  geometry)  is 
independent  of  the  spatial  y  variable,  and  thus  the  problem  can  be  worked  out  entirely  on  a  plane. 
The  scope  of  the  next  section  is  the  determination  of  the  EM  fields  within  region  2 — that  is,  the 
fields  in  the  virtual  LHM.  Once  the  fields  within  each  region  are  known,  pertinent  quantities  of 
interest  can  be  readily  evaluated:  the  spatial  electric  current  distribution  that  accounts  for  the 
electric  loss  mechanism,  the  polarization  current  induced  in  the  lossy  dielectric,  and  the  total 
radiated  field,  among  others. 

Following  a  standard  MAS  formulation,  two  equivalent  simulated  situations  are  set  up  to 
mimic  the  original  ones  in  regions  1  and  2.  In  region  1,  the  scattered  field  is  simulated  by  the  total 
field  generated  by  a  set  (set  I)  of  auxiliary  sources,  which  are  fictitious  y-directed  current 
filaments.  As  shown  in  Figure  4.4,  these  filaments  are  homogeneously  distributed  along  an 
auxiliary  surface  parallel  to  the  plane,  at  positions  given  by  p\i  =  xix-d1z  for  the  z-th  auxiliary 

source  of  set  I,  where  d\  is  the  distance  between  the  physical  interface  and  the  auxiliary  surface 
on  which  the  current  sources  of  set  I  lie.  These  filaments  carry  unknown  currents  and  are  treated 
as  sources  radiating  in  an  unbounded  free  space.  The  unknown  scattered  electric  field  at  the 
observation  point  in  region  1  is  described  by  the  superposition  of  the  radiated  EM  fields  of  the 
auxiliary  sources  of  set  I: 


Figure  4.4.  MAS  geometry  for  a  uniform  half-space  model. 
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(4.13) 


E°bJ  (x,z)  =  -  — aHr)(k  J(x  -  x  )2  +(z  +  d  )2 )  . 

4  V 

Accordingly,  the  corresponding  magnetic  field  components  are  expressed  as 


H°xbj(x,z)  =  jk0^~ 


z  +  d, 


4<s](x-x.)2  +  (z  +  d{)2 

a.H(2)  (ko  sjix  -  x.  )2  +  (z  +  d]  )2  ) 


(4.14) 


H°zbJ(x,z)  =  -jk0Yj 


T  4^j  (x-x.)2  +(z  +  dxf 


a.H(2)(k 


X.)2  +(z  +  dx)2) 


(4.15) 


Similarly,  in  the  dissipative  half-space  of  region  2,  as  shown  in  Figure  4.4,  the  field  is  simulated 
by  the  superposition  of  the  radiated  fields  of  a  set  (set  II)  of  bf  auxiliary  sources,  which  are  also  y- 
directed  current  filaments  and  are  once  again  homogeneously  distributed  along  an  auxiliary 
surface  parallel  to  the  plane.  The  z-th  auxiliary  source  of  set  II  is  located  at  p2i  =  X-X  +  d2z , 

where  d2  is  the  distance  between  the  physical  interface  and  the  auxiliary  surface  containing  the 
auxiliary  current  sources  of  set  II.  The  filaments  carry  unknown  currents  and  are  treated  as 
sources  radiating  in  an  unbounded  space  filled  with  the  LMH  medium  of  region  2.  The  EM  field 
within  region  2  is  then 


ELyHM(x,z)  =  ~^^y^biHl2\klJdcEcy  +(z-c?2)2  )  . 


The  corresponding  magnetic  field  components  in  the  LHM  can  be  expressed  as 


z  —  d  „ 


‘  4 ^/(x-x.)2  +(z-42)2 

a,H™  (k,  V(*-x.)2+(z-rf2)2) 


(4.16) 


(4.17) 


and 


■  4  ^/(x-x)2  +(z-42)2 

b,H™  (*,  7(x-xf)2+(z-42)2) 


(4.18) 
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where 


kx  =  k0J?™JJA™  (4.19) 

is  the  wave  number  in  the  LH  material  of  region  2  and  depends  on  the  relative  dielectric 
permittivity  and  magnetic  permeability.  The  relation  between  the  EM  fields  inside  and  outside  the 
LH  medium  is  dictated  by  the  boundary  conditions  (4.10).  More  specifically,  the  tangential 
components  of  the  electric  and  magnetic  intensities  must  be  continuous  across  the  boundary 
plane.  The  unknown  complex  currents  {«/}  and  i  =  1,  2,  3,  ...,  TV,  can  then  be  calculated  by 
imposing  the  boundary  conditions  (4.10)  on  2N  collocation  (matching)  points  on  the  interface. 
The  enforcement  of  the  boundary  conditions  at  z  =  0  leads  to  the  following  equations: 

H?  {kx  +  X  AP?  (K  V(xm-x,.)2+^2) 

(4.20) 

-  Y  BiHo2)  (kl  V(Xm“Xz)2+J 22  )  =  0 


for  m  =  1,  2,  3,  . . .,  M,  and 

H\2) (*o  +  ^  )  + 

xm  +h 

Y Ai  7  ^2 ,  =  (*» V^_x.)2+Ji2 >  (4-2 1 ) 

i  J(x  -x.)  +  d, 

-  /  ~^2^2  2gi2)(*.  V(^-^)2+d22)  =  0 

+j2 

for  m=  1,  2,  3,  ...,  M,  where  Ai  =  ai/Io  and  Bi  =  bi/Io  are  the  normalized  currents  of  the 

filamentary  fictitious  current  sources  of  sets  I  and  II  with  respect  to  the  current  of  the  primary 

source  /0.  The  auxiliary  sources  of  sets  I  and  II  and  the  collocation  points  are  taken  to  be  evenly 
spaced  on  their  respective  surfaces.  The  two  boundary  conditions  (4.10)  are  strictly  imposed  at 
M=  Nor,  alternatively,  at  M>  N selected  points  on  the  plane,  in  which  case  the  solution  of  (4.10) 
is  correct  in  a  least-squares  sense.  Usually  we  consider  the  former  situation,  in  which  the  number 
of  matching  points  is  identical  to  the  number  of  auxiliary  sources;  this  leads  to  a  linear  matrix 
equation  that  can  be  subsequently  solved  (via  standard  square-matrix  inversion  techniques)  for 
the  2N  unknown  quantities  At  and  which  correspond  to  the  unknown  currents.  Once  the  linear 
system  is  solved,  the  EM  fields  in  each  region  are  readily  computable  using  expressions  (4.13)- 
(4.19). 
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Actual 


Predicted  in  LHM 


IF  “  ^ 


Figure  4.5.  Electric  field  distribution  around  a  point  source  (left)  and  MAS-predicted  electric 
field  distribution  in  a  virtual  LHM  half-space  with  n  =  - 1  (right). 


4.6.  Numerical  results 

In  this  section  we  present  numerical  results  that  illustrate  the  capabilities  of  the  LHM 
approach  for  estimating  the  location  of  a  buried  object.  A  computational  code  has  been  developed 
using  the  formulation  of  the  previous  sections  and  assuming  equal  numbers  of  inner  and  outer 
auxiliary  sources.  Imposing  the  EM  boundary  conditions  at  matching  points  with  z  =  0  results  in 
a  square  matrix  system.  Both  the  auxiliary  sources  and  the  collocation  points  are  uniformly 
distributed  along  their  respective  surfaces.  In  our  subsequent  analysis  we  will  assume 
£^hm  =  jULHM  =  -1.  First,  numerical  tests  are  done  in  the  electromagnetic  wave  regime,  where  the 
target’s  size  is  comparable  to  the  relevant  wavelength.  Obviously  this  situation  is  very  different 
from  that  encountered  in  UXO  detection  and  discrimination  using  current  EMI  technologies. 
Therefore,  at  the  end  of  this  chapter  we  mention  several  steps  that  we  took  to  extend  the 
capabilities  of  the  LHM  approach  to  the  EMI  regime  and  discuss  its  practical  implementation  for 
realistic  applications. 


4.6.1.  Point-source  and  Gaussian-beam  illumination  results 

We  consider  point-source  and  Gaussian-beam  illumination  cases  to  study  whether  a 
planar  LH  medium  can  refocus  the  incident  beam  or  not.  We  first  investigate  the  point-source 
refocusing  problem.  A  point  source  placed  at  a  distance  h  =  30  cm  produces  a  field  on  the 
boundary  between  free  space  and  the  fictitious  LH  material  which  we  take  as  the  incident  field  on 
the  interface.  Using  the  boundary  conditions  we  can  determine  the  amplitudes  of  the  internal  and 
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external  auxiliary  sources  and  then  use  those  to  simulate  the  EM  fields  inside  the  virtual  LHM 
space.  Figure  4.5  shows  the  actual  field  from  a  point  source  (left)  and  the  field  predicted  using  the 
MAS  responding  amplitudes  (right).  Comparing  the  two  panels  we  can  see  that  the  predicted  field 
is  a  mirror  image  of  the  actual  source  field,  up  to  the  source  location.  The  distance  between  the 
interface  and  the  original  and  image  sources  is  virtually  the  same. 

Next  we  consider  a  Gaussian  beam  impinging  obliquely  on  the  boundary.  The  beam  was 
constructed  using  a  complex  radiation  source.  The  amplitude  of  the  incident  field  varies  spatially 

as  exp(-(x2  +  y2)/  w)  and  has  many  wave  vectors  associated  with  it.  The  wave  vectors  off  the 
beam  axis  point  away  from  the  axis  for  a  diverging  beam  and  towards  it  for  a  converging  beam.  It 
was  expected  that  a  LHM  with  negative  index  of  refraction  would  focus  the  beam;  i.  e. ,  it  would 
bend  the  wave  vectors  of  a  diverging  beam  back  towards  the  beam  axis.  The  electric  field 
intensity  produced  by  the  Gaussian  beam  source  in  free  space  is  shown  on  the  left  panel  of 
Figure  4.6.  The  highest  values  are  colored  red;  the  lowest  values  are  blue.  The  beam  is  clearly 
diverging  as  it  propagates  away  from  the  origin. 


Figure  4.6.  Electric  field  distribution  from  a  Gaussian  beam  (left)  and  MAS-predicted  electric 
field  distribution  in  a  virtual  LHM  half-space  with  n  =  - 1  (right). 


73 


Figure  4.7.  Four  Gaussian  beam  sources  placed  in  a  virtual  LHM  space  with  n  =  - 1.  Left:  electric 
field-intensity  distribution.  Right:  phase  distribution. 


The  right-hand  panel  of  Figure  4.6  shows  the  reconstructed  beam  inside  a  LHM  half¬ 
space  with  n  =  - 1  after  solving  the  full  EM  problem.  The  beam  source  in  the  simulation  was 
placed  at  a  distance  of  50  cm,  and  the  wavelength  is  3  cm;  the  interface  boundary  ranged  between 
-5  m  and  5  m.  The  predicted  result  clearly  shows  that  the  LH  planar  medium  turns  the  diverging 
wave  vectors  toward  the  beam  focus  point  and,  hence,  acts  as  a  lens  that  focuses  the  beam.  As 
mentioned  in  the  previous  paragraph,  the  distance  between  the  interface  boundary  and  the  original 
source  coincides  with  that  between  the  boundary  and  the  reconstructed  source.  This  means  that 
the  LHM  half-space  ideally  reproduces  the  actual  beam  in  the  high-frequency  EM  regime. 


4.6.2.  LHM  applied  to  multiple  sources  and  objects 

In  this  section  the  LHM  approach  is  applied  to  scenarios  involving  multiple  objects.  We 
first  consider  the  reconstruction  of  four  Gaussian  beam  sources  enclosed  within  a  virtual  LHM 
with  n  =  - 1.  The  MAS  technique  is  used  to  solve  the  EM  boundary  value  problem  and  predict  the 
EM  field  in  the  LHM  space.  The  resulting  field-intensity  and  phase  distributions  are  shown  in 
Figure  4.7.  The  field-intensity  distribution  clearly  shows  that  the  LHM  refocuses  the  EM  field 
from  the  individual  targets  without  losing  information  on  the  individual  scatterers.  In  addition,  the 
phase  distribution  illustrates  that  the  LHM  technique  requires  very  accurate  matching  conditions 
on  the  boundary.  This  indicates  that  the  technique  could  potentially  face  problems  in  the  detection 
and  discrimination  of  subsurface  objects  in  noisy  environments  (as  is  the  case  in  any  realistic 
UXO-related  application). 
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Figure  4.8.  Actual  (in  free  space,  at  left)  and  reconstructed  (in  a  LHM  with  n  =  - 1,  at  right) 
electric  field-intensity  (top  row)  and  phase  (bottom  row)  distributions  for  three  point 
sources. 


Figure  4.9.  Actual  (in  free  space,  at  left)  and  reconstructed  (in  a  LHM  with  n=—  1,  at  right) 
electric  field  intensity  (top  row)  and  phase  (bottom  row)  distributions  for  a 
cylindrical  scatterer. 

To  further  illustrate  the  capability  of  the  virtual  LHM  as  a  way  of  detecting  and 
separating  radiating  objects  we  consider  three  point  sources  placed  in  free  space  and  separated  by 
one  wavelength.  The  LHM  reconstruction  algorithm  is  implemented  as  before:  the  EM  fields 
from  the  sources  are  calculated  on  the  boundary  between  free  space  and  the  virtual  LHM.  No 
other  information  is  given.  By  solving  a  linear  least-squares  problem  we  can  reconstruct  the  EM 
fields  inside  the  virtual  LHM  space.  The  actual  (left)  and  reconstructed  (right)  distributions  of 
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Figure  4.10.  Electric  field  intensity  distribution  for  three  objects:  actual  (in  free  space,  at  left) 
and  reconstructed  (in  a  LHM  with  n  =  -1,  at  right). 


Figure  4.11.  Magnetic  field  intensity  (7/x)  distribution  for  two  objects:  actual  (in  free  space,  at 
left)  and  reconstructed  (in  a  LHM  with  n  =  - 1,  at  right). 
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electric  field  intensity  (top  row)  and  phase  (bottom  row)  are  depicted  in  Figure  4.8.  The  results 
clearly  show  that  the  proposed  technique  is  capable  of  inferring  the  locations  of  multiple  objects 
if  the  scatterers  are  separated  by  at  least  one  wavelength.  The  analysis  has  been  done  in  the  wave 
regime  while  keeping  in  mind  that  it  should  eventually  be  reimplemented  in  the  low-frequency 
EMI  regime  if  it  is  to  be  applied  to  the  UXO  detection  and  discrimination  problem.  Next  we  study 
the  EM  scattering  from  an  extended  object,  a  cylinder.  The  results  appear  in  Figure  4.9. 

In  our  next  example  we  reconstruct  the  EM  field  due  to  three  objects  separated  by  more 
than  one  wavelength  and  illuminated  by  a  Gaussian  beam.  The  results,  given  on  Figure  4.10, 
depict  the  actual  electric  field  intensity  (left)  and  the  reconstructed  intensity  in  the  virtual  LHM 
(right).  The  LHM  reconstruction  algorithm  is  able  to  infer  the  targets’  locations  and  predicts  very 
accurately  the  EM  field  distribution  between  the  free-space/LHM  boundary  and  the  objects.  On 
the  other  hand,  the  reconstructed  fields  behind  the  objects  diverge  from  the  actual  field;  this  is 
because  the  scattered  field  at  the  free-space/LHM  interface  does  not  contain  information  about 
the  EM  field  due  to  the  objects’  invisible  parts. 

Finally,  we  apply  the  LHM  reconstruction  technique  to  a  cylindrical  object  immersed  in 
another  cylinder.  The  targets  are  illuminated  by  a  Gaussian  beam  source;  the  immersed  object  has 
a  size  equal  to  one  wavelength.  The  scattered  EM  field  from  the  targets  is  computed  on  the 
surface  of  the  LHM  half-space  and  then  reconstructed  inside  the  medium.  The  actual  and 
reconstructed  Hx  and  Hy  magnetic  fields  are  depicted  on  Figures  4.1 1  and  4.12  respectively.  The 
results  show  that  the  algorithm  is  able  to  detect  and  correctly  identify  the  scatterers’  geometries 
and  positions  using  different  components  of  the  scattered  magnetic  field. 


Predicted  in  LHM 


Figure  4.12.  Magnetic  field  ( Hy )  intensity  distribution  for  two  objects:  actual  (in  free  space,  at 
left)  and  reconstructed  (in  a  LHM  with  n  =  - 1,  at  right). 
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Figure  4.13.  Reconstructed  electric  field-intensity  (left)  and  phase  (right)  distributions  for  three 
point  sources  in  a  LHM  with  n  =  - 1  when  the  wavelength  is  1000  times  greater 
than  the  distances  between  the  objects. 


4.6.3.  An  attempt  to  extend  the  LHM  approach  to  the  EMI  regime 

In  the  EMI  (low-frequency)  regime  the  relevant  wavelengths  are  much  larger  than  the 
sizes  of  the  scatterers,  and  as  a  consequence  there  are  no  waves,  but  rather  diffusion  and 
potential  fields.  In  order  to  illustrate  the  performance  of  the  LHM  approach  in  the  low-frequency 
case  we  revisit  the  problem  reported  in  Figure  4.8  but  make  the  wavelength  be  1000  times  greater 
than  the  distance  between  the  objects.  The  results  appear  in  Figure  4.13.  Now  the  three  objects 
seem  to  behave  as  a  single  target,  and  it  becomes  very  difficult  to  identify  their  locations 
correctly. 

This  suggests  that  we  cannot  directly  apply  LHM  wave  phenomena  to  the  EMI  case,  as 
there  is  no  field  propagation  per  se.  However,  we  thought  that  by  using  a  spatial  distribution  of 
EMI  measured  data  we  could  generate  a  synthetic/theoretical  Gaussian  beam.  As  is  well 

documented  in  the  optics  literature,  a  Gaussian  (paraxial)  beam  is  a  plane  wave  e~jkz  (with  wave 
number  k  =  2n!X  and  wavelength  A)  modulated  by  a  complex  envelope  A( r)  that  is  a  slowly 
varying  function  of  position.  The  complex  amplitude  is  thus 


U(r)  =  A(r)  e  jkz .  (10) 

The  envelope  is  assumed  to  be  approximately  constant  within  a  neighborhood  of  size  A,  so  the 
beam  behaves  locally  like  a  plane  wave.  The  spatial  distribution  of  the  envelope  A( r)  determines 

the  beam  shape  (Gaussian,  rectangular,  sine  or  cosine)  and  the  factor  e-jkz  describes  wave 
propagation.  In  EMI  we  cannot  incorporate  this  latter  term,  but  we  can  still  try  to  match  the 
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Figure  4.14. 


Problem  geometry  (left)  and  reconstructed  magnetic-field  intensity  distribution 
(right)  for  a  dipole. 


spatial  distribution  of  the  envelope  A( r)  to  monostatic  or  bistatic  data  D( r)  measured  on  a  set  of 
grid  points;  this  can  be  either  (l)the  magnetic  field  Hz  measured  by  the  GEM  sensor,  (2)  the 
induced  voltage  in  the  EM-63  sensor,  or  (3)  the  amplitude  of  the  total  scattered  magnetic  field  as 
measured  by  the  new  GEM-3 D  sensor. 

To  demonstrate  the  applicability  of  the  proposed  technique  to  the  EMI  problem  we  focus 
on  a  simple  example.  Above  a  LHM  slab  we  place  a  static  magnetic  dipole  oriented  along  the  z- 
axis  (see  Figure  4.14).  Using  the  Biot-Savart  law  we  can  compute  the  spatial  distribution  of  the 
Hz(y)  magnetic  field  on  the  top  boundary  of  the  LHM  slab  and  take  it  to  be  the  shape  of  the 
envelope:  A( r)  =  Hz( r).  We  then  solve  the  full  EM  scattering  problem  for  normal  incidence  of  the 
Gaussian  beam  on  the  slab.  The  reconstructed  near-field  distributions  inside  and  outside  the  LHM 
are  shown  in  Figure  4.14.  It  is  worth  noting  that  the  distance  from  the  top  boundary  of  the  LHM 
to  the  first  focus  (white  arrow)  equals  the  distance  from  the  first  focus  to  the  dipole  (black  arrow). 
This  illustrates  that  LHM  techniques  could  potentially  be  used  to  locate  buried  objects.  However, 
in  this  test  we  chose  a  “synthetic”  wavelength  A  comparable  to  the  distance  between  the  matching 
boundary  and  the  dipole.  Further  studies  demonstrated  that  the  predicted  dipole  location  depends 
on  A,  and  that  in  general  the  reconstructed  locations  for  low-frequency  cases  are  very  sensitive  to 
the  actual  locations  and  to  the  spatial  range  and  density  of  the  measured  data.  The  least-squares 
problems  (4.20)  and  (4.21)  are  also  sensitive  to  the  noise  level.  All  of  these  facts  make  the  LHM 
approach  less  attractive  for  UXO  detection  and  discrimination. 
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4.7.  Conclusion 


In  this  chapter  we  investigated  the  LHM  approach  as  a  potential  candidate  for  estimating 
the  location  of  a  buried  object  without  having  to  solve  a  traditional  ill-posed  inverse  scattering 
problem.  We  discussed  in  detail  the  theoretical  basis  and  the  numerical  realization  of  the  method 
in  the  high-frequency  EM-wave  regime  and  then  used  it  to  reconstruct  the  field  and  infer  the 
location  of  a  point  source,  a  Gaussian  beam,  and  different  combinations  of  multiple  objects.  The 
results  demonstrated  the  refocusing  property  of  a  planar  LHM  half-space  (/.£.,  the  fact  that  it 
sends  signals  backward  to  the  origin)  and  its  ability  to  estimate  the  location  and  geometry  of  an 
object.  These  results  suggest  a  number  of  interesting  applications  of  the  LHM  approach.  Lor 
example,  it  could  be  used  as  a  new  tool  for  solving  inverse  scattering  problems  in  the  microwave 
and  optical  frequency  regimes.  Another  potential  application  could  be  image  reconstruction  with 
sub-wavelength  resolution. 

We  then  applied  the  proposed  technique  to  low-frequency  cases  in  order  to  assess  its 
potential  usefulness  in  the  context  of  the  UXO  problem.  Our  studies  show  that  in  the  low- 
frequency  regime  the  LHM  refocusing  approach  is  impractical  because  the  necessary  wavelength 
is  much  larger  than  the  size  and  depth  of  any  potentially  interesting  target.  We  attempted  to 
construct  a  Gaussian  beam  by  combining  a  monostatic  EMI  spatial  distribution  with  an 
empirically  assigned  wave-phase  part.  However,  our  studies  reveal  that  the  approach  is  very 
sensitive  to  the  object’s  orientation  and  to  the  wave  number.  Thus  the  refocused  EM  fields  have 
maxima  whose  locations  depend  quite  strongly  on  the  phase  distribution  on  the  boundary.  In 
addition  we  found  that  the  method  is  very  sensitive  to  noise.  The  fact  that  real  UXO  detection  and 
discrimination  problems  always  involve  very  noisy  data  leads  us  to  conclude  that  the  LHM 
approach  is  more  suitable  for  the  EM  wave  regime  than  it  could  be  for  UXO  discrimination. 
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Chapter  5 

Discussion  and  Conclusions 


5.1.  Objectives 

The  objective  of  the  proposed  work  was  to  develop  new  non-traditional  physics-based 
inverse  approaches  that  could  let  us  determine  the  location  and  orientation  of  buried  objects 
starting  from  actual  EMI  data  but  without  solving  an  ill-posed  inverse  scattering  problem.  The 
specific  technical  objectives  were  the  following: 

1.  Develop  an  approach  combining  NSMS  and  the  physics-based  pole  distribution  and 
prove  that  this  combination  can  be  used  to  determine  the  location  and  orientation  of  a 
buried  object  starting  from  measured  spatially  distributed  EMI  data  without  having  to 
resort  to  the  traditional  inverse-scattering  procedure. 

2.  Determine  if  the  directional  vectors  corresponding  to  the  energy  gradient  of  the  scattered 
magnetic  field  always  cross  at  a  single  point,  and  whether  or  not  this  point  is  always 
inside  the  scatterer.  Use  NSMS  to  perform  upward  continuation  on  the  measured  data  in  a 
defined  space,  compute  the  magnetic  field  gradient  at  each  point  in  the  space,  determine 
the  gradient’s  directional  vectors  toward  the  ground,  find  the  point  where  these  vectors 
cross,  and  relate  the  crossing  point  to  the  buried  object’s  location. 

3.  Determine  the  feasibility  of  exploiting  the  physical  properties  of  new  left-handed  media 
(LHM),  in  particular  the  possibility  of  making  perfect  lenses  out  of  them,  to  determine 
the  location  and  orientation  of  an  object.  Construct  a  paraxial  (Gaussian)  beam  with  the 
spatial  part  of  an  actual  set  of  EMI  measurements,  define  a  left-handed-medium  slab,  and 
refocus  the  beam  inside  the  slab,  thus  determining  the  object’s  location. 

4.  Evaluate  and  document  the  robustness  of  the  proposed  approaches  for  most  realistic 
scenarios,  in  which  an  object’s  response  is  contaminated  by  complex,  heterogeneous 
noise  associated  with  the  positioning  accuracy  of  the  sensor,  uncertainties  in  orientation, 
and  the  possibility  of  response  from  cluttered  background  involving  two  or  more  objects. 

Under  this  project  we  developed  three  methods:  (1)  HAP,  (2)  Pole-series  expansions  and 
(3)  LHM-based  approaches.  We  assessed  the  ability  of  each  technique  to  predict  the  location  and 
orientation  of  a  buried  object  by  comparing  inversion  results  to  the  ground  truth.  We  found  that 
all  three  techniques  have  the  potential  to  estimate  position  and  orientation  at  different  frequencies. 
However,  when  objects’  responses  are  contaminated  with  additive  noise,  the  performance  of  the 
pole-series-expansion  and  LHM  methods  degrade  significantly  compared  to  that  of  the  HAP 
technique.  In  addition,  the  computational  time  required  by  the  PSE  and  LHM  approaches  and 
similar  practical  considerations  prevented  us  from  further  testing  and  developing  those  two 
techniques,  leaving  the  HAP  as  our  weapon  of  choice  for  determining  buried  objects’  locations 
and  orientations. 
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5.2.  HAP  method 


The  HAP  method  is  a  new  physics-based  analytical/numerical  technique  that  uses  global 
quantities — the  magnetic  field  H,  the  vector  potential  A,  and  the  scalar  potential  y/- — to  estimate 
the  location,  orientation,  and  magnetic  polarizability  of  a  buried  object.  The  method  should  be 
valid  for  both  land-based  and  underwater  UXO  surveys.  To  find  A  and  y/  from  measured  H  (in 
monostatic  or  multistatic  cases)  we  developed  a  two-dimensional  version  of  the  normalized 
surface  magnetic  source  (NSMS)  model,  in  which  fictitious  sources  are  distributed  on  an 
auxiliary  surface  just  below  the  measurement  grid  and  their  responding  amplitudes  are 
determined  using  a  collocation  algorithm.  We  tested  the  HAP  technique  in  the  time  and  frequency 
domains,  using  both  synthetic  and  measured  data,  in  noise-free  and  noisy  environments,  and  with 
point  dipoles,  spheres,  and  actual  UXO  as  scatterers. 

We  carried  out  blind-test  analyses  using  real  UXO  interrogated  by  the  MPV-TD  and 
GEM-3D+  sensors.  We  first  estimated  the  positions  and  orientations  of  the  objects  via  HAP,  and 
then  to  perform  classification  we  used  a  three-dimensional  NSMS  method  with  dipole  responding 
sources  distributed  on  a  spheroid  surrounding  the  object:  for  each  unknown  target  we  determined 
the  NSMS  amplitudes,  integrated  those  to  find  a  total  time/frequency-dependent  NSMS  curve, 
and  used  the  curves  as  discriminators  by  comparing  them  to  pre-computed  samples  belonging  to  a 
library  of  known  UXO  types.  The  method  was  able  to  identify  with  reasonable  accuracy  the 
positions  and  orientations  of  the  targets  and  had  100%  success  classifying  them  into  UXO  and 
clutter  and  identifying  the  UXO  type. 

We  also  adapted  the  HAP  method  to  the  EM63  and  NRL  EMI  array  TD  sensors.  Our 
results,  reported  in  Chapter  2,  clearly  demonstrate  the  superior  performance  (and  the  robustness 
with  respect  to  sensor  noise)  of  the  HAP  technique  as  a  method  for  estimating  the  locations  and 
orientations  of  single  or  aggregate  buried  targets. 


5.3.  Outlook 

The  work  presented  here  has  demonstrated  that  the  HAP  approach  has  the  potential  to 
make  UXO  discrimination  more  reliable  and  faster  by  providing  a  trustworthy  means  to 
determine  location  and  orientation.  However,  a  substantial  research  and  development  effort  must 
be  undertaken  before  the  technique  can  be  applied  to  real-world  UXO  problems.  This  effort 
should  address,  among  others,  the  following  issues: 

•  Comparison  to  nonlinear  inverse-scattering  models:  Nonlinear  inverse-scattering 
models  are  the  most  commonly  used  techniques  for  determining  the  position  of  a  buried 
object.  In  this  study,  we  did  not  investigate  if  the  HAP  approach  represents  a  significant 
improvement  over  traditional  nonlinear  inverse-scattering  methods.  We  have  found  HAP 
to  be  fast,  accurate,  and  not  plagued  by  local  minima,  which  suggests  that  the  technique 
can  be  used  in  combination  with  physically  complete  models — such  as  NSMS  and 
GSEA — to  perform  UXO  discrimination  in  real  or  near  real  time.  We  would  expect  that, 
once  an  object’s  position  and  orientation  are  accurately  known,  a  UXO  discrimination 
approach  based  on  physically  complete  forward  models  would  be  much  faster  and  better 
than  the  traditional  nonlinear  inverse  models  because  the  target’s  intrinsic  parameters 
could  be  determined  much  more  exactly.  Our  studies  show  that  the  HAP  model  is  robust 
with  respect  to  noise.  Now  that  we  have  established  a  method  for  recovering  an  object’s 
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location  and  orientation,  we  should  be  in  a  position  to  compare  our  algorithm  to  more 
traditional  approaches  with  regard  to  stability  and  generality. 

•  EMI  data  density  and  sensitivity  analysis:  For  application  to  practical  UXO 
discrimination  problems,  it  is  critical  to  determine  the  data  fidelity  required  for  HAP 
techniques.  In  order  to  estimate  the  2D  NSMS  responding  amplitudes  and  the 
polarizability  tensor  of  the  induced  magnetic  dipole  we  need  to  determine  a  minimum 
number  of  data  points  (or  a  minimum  data  density)  as  a  function  of  signal-to-noise  ratio, 
spatial  data  coverage,  positional  error,  and  other  survey/instrument  characteristics. 

•  Multiple  Targets:  The  data  considered  in  this  study  were  from  single,  isolated  targets. 
However,  separating  and  performing  discrimination  on  overlapping  signals  from  multiple 
objects  using  the  HAP  or  any  other  existing  technique  is  still  an  open  question.  To  extend 
HAP  into  an  A-target  locator  technique,  one  that  will  reliably  say  how  many  targets  are 
present  in  the  sensor’s  field  of  view,  where  they  are,  and  what  tilt  they  have,  all  of  this 
without  resorting  to  computationally  taxing  optimizations,  is  a  highly  desirable  and 
important  step  that  the  UXO  community  must  take. 

•  Geologic  background:  Signals  coming  from  background  geology  always  pose  additional 
problems  during  UXO  discrimination.  We  need  to  investigate  how  the  presence  of 
magnetically  susceptible  soil  affects  the  performance  of  HAP;  i.  e. ,  we  need  to  incorporate 
the  soil’s  response  into  the  HAP  model. 

•  Underwater  UXO  discrimination:  We  have  demonstrated  mathematically  that  the  HAP 
approach  can  in  principle  make  valid  predictions  if  the  target  is  in  an  underwater 
environment.  At  this  point,  however,  we  have  not  investigated  if  this  is  true  when  actual 
EMI  data  are  considered.  Our  success  with  land-based  predictions  is  a  good  first  step  in 
that  direction. 
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